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ABSTRACT 


The response of a well defined impervious surface to a specified 
input of rainfall was studied by means of a laboratory catchment. Under 
natural conditions these responses are affected by many random factors 
such as evaporation and infiltration losses. Under laboratory conditions 
such factors are eliminated or negligible. Thus the laboratory experi- 
ments enable investigation of deterministic processes in hydrology. 

A mathematical model of overland flow has been developed © i 
using the state variable approach to system analysis. The model was 
tested by simulating the experimental hydrographs of surface runoff from 
the laboratory catchment for uniform, and time and space variable rain- 
fall inputs. A mathematical relationship describing the raindrop impact 
effect was developed by the method of dimensional analysis and subjected 
to empirical evaluation. 

The experimental hydrographs were also simulated by two existing 
mathematical models of overland flow: the instantaneous unit hydrograph 
and the kinematic wave model. The purpose of the investigation was to 
study nonlinearity of the rainfall-runoff relationship, to examine the 
theoretical assumptions of the above mathematical models, and to provide 
a basis for comparison of the results obtained by the state variable 
model. 

The results of experiments indicated that the rainfall-runoff 
process is highly nonlinear, and thus the linear system analysis methods 


have limited application. The kinematic wave model was found to provide 
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fairly consistent results but has a tendency to over-estimate the time 

to peak and the peak discharge. Results obtained with the state variable 
model compared favourably with those obtained experimentally by means of 
the laboratory catchment. It is believed that the output of this 
mathematical model may be helpful in urban runoff calculations, and that 
the concept of state and state variables should find a wide application 


in general hydrologic simulation. 
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fairly consistent results but has a tendency to over-estimate the time 

to peak and the peak discharge. Results obtained with the state variable 
model compared favourably with those obtained experimentally by means of 
the laboratory catchment. It is believed that the output of this 
mathematical model may be helpful in urban runoff calculations, and that 
the concept of state and state variables should find a wide application 


in general hydrologic simulation. 
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CHAPTER I 


INTRODUCTION 


1-1 History 


Natural sciences are observational sciences. The first flood 
observations which were made, according to biblical accounts about 
2957 B.C. on the Euphrates, may thus properly indicate the beginning of 
hydrology. However, it was not until the 20th century A.D. that hydrology 
developed into a truly scientific discipline. 

In antiquity it was considered that stream flow was supplied from 
huge subterranean reservoirs and that water somehow returned from the 
oceans to the reservoirs through subterranean channels. No direct link 
was considered to exist between rainfall and runoff. 

The first recorded attempts at measuring the adequacy of precipi- 
tation to provide runoff were made in the 17th century by Perrault (1674) 
and Mariotte (1686), who measured the flow and precipitation of the 
Seine basin in France and concluded that precipitation falling in this 
basin was sufficient to supply the stream flow. 

From then onwards, in the 17th and 18th centuries, many obser- 
vations were made of hydrological phenomena and the period was marked 
by the development of important measuring equipment, such as the current 
meter for the measurement of velocity in streams. 

From the beginning of the 19th century to about 1930, hydrology 
went through the period of coefficients. This was, in a way, an 
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empirical approach, in which "universal propositions were made about 
causal connections" (Toebes, 1968). Propositions were usually expressed 
in a simple equation in which one or more coefficients expressed the 
universality and had to be selected by judgement based on experience. 

The excessive use of simple formulae, the results of which were 
not always found satisfactory, forced engineers to increase hydrological 
investigations to find more rational methods, based on a reasoning pro- 
cess rather than on judgement. 

From about 1930 hydrology went through a period of great advance. 
The hydrological cycle in all its aspects became the central concept of 
the science, and rather than the restricted view of hydrology held before 
as a subject largely dealing with stream flow, it became a science which 
dealt with water in all its aspects. Definitions of hydrology such as 
"Hydrology is an earth science which is concerned with the laws governing 
the time and space distribution, circulation, and properties of the 
waters of the earth" were made and the sphere of hydrological interest 
has become interrelated in many aspects with other sciences, ranging 


from meteorology to chemistry and geology. 


1-2 Rainfall - Runoff Relationship 


Engineers are often concerned with only part of the hydrologic 
cycle which manifests itself as river runoff. Because rainfall is a 
prime cause of runoff in many parts of the earth, understanding of the 
rainfall-runoff relationship is essential for engineering analysis. The 
present emphasis on the development of water resources accelerates the 
need to develop reliable theories for the solution of problems such as 


the prediction of the future runoffs from watersheds altered by man-made 


once s: seeks to vers B dawoxn sob rate ar soot ne 
‘Re Iqson02. hake cs ani ‘amabad pdaaqee ast He ab efate & 
sented bad vanloxbys to wexby besalxI pos, eiy stig ts 


Aoki eonat a b) amkbed’ db vical densi dgtwg 


changes, the design of a reservoir capacity for irrigation or water- 
power development, and the prediction of flood flows in rivers. 
Although a large number of theories and methods now exists, there is a 
general lack of confidence concerning their ability to produce reliable 
solutions to the above type of engineering problem. 

River runoff phenomena, besides their relation to rainfall as 
their sources, are the result of a highly complex combination of pro- 
cesses taking part at the surface and in layers of the ground. These 
processes in general are non-stationary because of temporal changes in 
the parameters governing them. Some of these parameters, such as soil 
moisture potential, soil permeability, and land-surface roughness, are 
variable over short periods of time. 

Like all physical phenomena on the earth the processes which 
transform rainfall into runoff conform to che laws of physics and in 
this sense are completely regular. However, the extraordinary variety 
of conditions under which these processes occur, the non-uniformity of 
their spatial and temporal development, and finally the random character 
of fluctuations and combinations of different factors cause the analysis 


of the rainfall-runoff relationship to be a difficult task. 


1-2-1 Traditional Analysis 


Classical hydrology has attempted to solve the relation between 
rainfall and runoff on a given catchment with records through a study of 
the dependence of the volume of effective rainfall in an individual 
event on the volume and duration of actual rainfall, the time of year, 
etc., usually expressed in co-axial graphical correlation. For attenu- 


ation the unit hydrograph method was used. For ungauged catchments the 
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traditional approach has been to relate the unit hydrograph to the 
catchment characteristics, such as the catchment area and shape, the 
length and slope of the main stream, etc. 

The unit hydrograph concept was proposed by Sherman (1932) as a 
means of transforming effective rainfall (i.e. rainfall from which all 
losses due to evaporation, deep seepage, etc., have been deducted) into 
a runoff hydrograph. The unit hydrograph of a drainage basin is defined 
as a hydrograph of direct runoff (i.e. that part of runoff which enters 
the stream promptly after rainfall) resulting from a unit depth of 
effective rainfall generated uniformly over the basin area at a uniform 
rate during a specified period of time. It is assumed that the direct 
runoff is derived from the effective rainfall by a linear operation. 
Hence, if the unit hydrograph for the drainage basin is known, the 
hydrograph due to any effective rainfall distribution can be derived by 
applying principles of proportionality and superposition of causes and 
effects, which are valid for linear processes. 

Derivation of the unit hydrograph from an observed hydrograph 
requires estimates of the depth and duration of the effective rainfall. 
Traditionally the observed total runoff hydrograph has been separated 
into the direct runoff and base flow. Parameters of the effective 
rainfall are then determined from the amount of the direct runoff. 
However, because the total runoff hydrograph consists of a continuum of 
components (surface runoff, subsurface runoff, and groundwater runoff) 
it is virtually impossible to distinguish them and thus the separation 
of the hydrograph becomes extremely arbitrary. This arbitrariness of 


the base flow separation necessarily destined further rainfall-runoff 
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analysis, and the relationships thus derived to be somewhat haphazard. 


1-2-2 New Approaches 


The strength of hydrology lies in the excellence of the basic and 
derived sciences on which it is founded and the capability with which 
hydrologists can apply findings of these sciences to solve hydrological 
problems. Both the rapid development of related sciences such as mathe- 
matics, engineering and computing sciences, and the concentrated effort 
during the International Hydrological Decade have prompted introduction 
of new techniques into the rainfall-runoff analysis. The majority of 
these techniques can be categorized into one of the following general 
classifications: (1) system approach to the rainfall-runoff relation- 
ship, and (2) laboratory catchment studies of the rainfall-runoff 


relationship. 


1-2-2-1 System Approach 
The difficulty of applying known physical laws to hydrological 


phenomena lies in the ignorance of boundary and input conditions (Nash, 
1969). Modern hydrology has attempted to circumvent these difficulties 
by a system approach. In the system approach, instead of studying the 
laws governing the operation and applying these to the boundary conditions, 
the operation itself, as manifest in previously recorded events, is 
studied. 

The concept of a system, borrowed from electrical engineering, 
found widespread use by hydrologists. A system can be defined (Naslin, 
1965) to be an assembly of physical components which experience mutual 


actions and reactions under the influence of their own behaviour and the 
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external excitations acting upon the system as a whole. Essentially, 

a system consists of three parts (Vemuri, 1970): (1) the input signal, 
(2) the output signal, and (3) the plant. The plant of a system contains 
a mathematical characterization of the process that relates the inputs 
and outputs. 

Problems solved by system engineering can be broadly classified 
as direct and inverse problems. Direct problems are characterized by 
a complete specification of a system and the task is to study or predict 
the response of the system to any specified input. The inverse problems 
are much more complex. Here the response to a particular input or inputs 
is known, but either the equations describing the process are unknown, 
or the inputs themselves are unknown. Depending upon the specific cir- 
cumstance, an inverse problem may present itself as (1) design problen, 
(2) control problem, and (3) identification or modeling problem. 

The task of analysing the rainfall-runoff relationship and 
simulating the runoff process by a mathematical model is a typical 
identification problem of systems analysis. In general, an identification 
problem is characterized by a specification of a finite set of observa- 
tions on the input and output signals of a system and the goal is to 
find a mathematical characterization of the system. There are two 
principal lines of attack in attempting to characterize a system (Vemuri, 
1970): 

(1) The dynamical (or differential) equation approach. The 
unknown system is assumed to be characterized by a known set of dynamic 
(often differential) equations of a given order but with unknown 
coefficients. State variable and transfer function methods belong to 


this class. 
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(2) The kernel function approach. No assumption is made about 
the unknown system which is characterized by an analytic function over a 
function space and represented by the associated power series expansion 
(Balakrishnan, 1963). Identification of the inputs response essentially 
belongs to this class. 

A basic advantage of the system approach lies in the possibility 
for standardization of analytical and computational techniques. Once a 
handful of techniques are mastered, it becomes only necessary to extend 
the arguments to be applicable to any physical system (Vemuri, 1970). 
Furthermore, the system representation helps to identify analogous 
situations which may otherwise escape attention. General as it may 
seem, the system approach is not a substitute for the lack of under- 
standing of the governing laws of the studied process. There are two 
aspects of the system approach which must be realized. 

Firstly, the form of the equation, or the model of a system, 
cannot be found through analysis of the input and the output (black box 
approach), but must be assumed. The system analysis says very little 
about the choice of the model, or about how a specific system should be 
characterized. Rather, it demonstrates the application of mathematical 
tools after the system has been properly represented by a mathematical 
model. 

Secondly, the identification problem may lead to a multitude of 
mathematically acceptable solutions out of which one physically acceptable 
solution, if any such solution exists, has to be selected. However, it 
is almost practically impossible to verify the validity of an identified 


model. Two main reasons for this are: (1) in an identification problem 
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the input is not generally under the control of the observer, and (2) the 
period of observation is always limited and therefore the validity of 

the model is correspondingly limited to the range of variation of input 
signals which appeared within the observation interval. The wording 
"identification of a system from input and output measurements" is 
therefore not strictly correct; the best one can do is to make an esti- 
mate, from the measurements, of the best "representative" in a given 


class (Balakrishnan, 1967). 


1-2-2-2 Laboratory Catchment Studies 


The conventional approach to study the rainfall-runoff relation- 
ship of a watershed uses historical data to fit a mathematical model for 
simulation of watershed hydrologic behaviour. Historical data, however, 
are collected from natural phenomena that can be measured and observed 
only once and then will not occur again. Furthermore, the available 
data usually cover too short a period and contain too limited a range of 
events to provide a thorough test of mathematical models of the runoff 
process. For these reasons there has been growing interest in the studies 
of the rainfall-runoff relationship under controlled conditions, by means 
of laboratory catchments. 

Amorocho (1965) defined a laboratory catchment as a physical 
system of comparatively small size, in which the input and the output 
occur in the form of flowing water, and appropriate controls and measure- 
ment devices are provided to regulate the input as desired and to measure 
the output with adequate accuracy. The majority of laboratory catchments 
can be categorized into one of the following general classifications: 


(1) model catchments, (2) hydromechanic prototypes, and (3) prediction 
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analysis prototypes. 

The first type of laboratory catchment is used in studies of the 
rainfall-runoff relationship, which involve construction of a scale 
model of an actual watershed. The aim of such studies is to predict the 
prototype behaviour from observations made on the model, according to 
Newton's law of similitude. Several studies with model catchments of 
this type have been reported (Grace, 1965, 1967; Chery, 1965, 1966, 
1967), however, the complexity of the hydrologic processes poses very 
difficult problems with regard to the establishment of dynamic simili- 
tude. Not only does it appear unlikely that even the most significant 
requirements of similitude may be satisfied simultaneously, but it is 
a hard task to characterize rigorously most watershed parameters. 

The second of the above classifications includes laboratory 
catchments which are used for experimental investigation of flow 
mechanisms and are part of the field of study sometimes called "watershed 
hydraulics" (Amorocho, 1965; Chow, 1967). With the use of this type of 
a laboratory catchment as a research tool, a number of studies on the 
rainfall-runoff relationship have been made, including studies of the 
time distribution of the runoff process (Harbaugh, 1966), the conceptual 
watershed roughness (Harbaugh and Chow, 1967), and the effect of the 
storm movement (Yen and Chow, 1968). There are numerous other problems 
that can be studied systematically using this experimental approach. The 
experimentation may include changes in physical variables of the catchment, 
such as the shape, roughness, slope, size, and channel pattern. The study 
may not be limited to the surface runoff only. The laboratory catchment 


may be built with perforated bottom and provided with a layer of porous 
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material. Thus, subsurface flow can be developed, tested, and analyzed. 
Information obtained from such experiments is otherwise unobtainable in 
the case of natural hydrological data. 

The type of a laboratory catchment depicted in the third classi- 
fication serves the purpose of providing laboratory data for the tests 
of mathematical models of the rainfall-runoff relationship built on the 
basis of the system approach, i.e., to provide data necessary for 
identification of the system. When this type of a laboratory catchment 
is used, no attempt is made to simulate the detailed inner behaviour of 
any specific watershed. The only requirement is that the laboratory 
catchment should have significant nonlinearities of the same kind 
presumed to exist in natural basins, and the prime purpose is to gain 
knowledge on methods for establishing the input-output relationship 
applicable to the general class of nonlinear systems of which watersheds 
are a part (Amorocho, 1965). Most experimental investigations dealing 
with this type of a laboratory catchment have been done by Amorocho and 


Orlob (1961), Amorocho and Hart (1965) and Amorocho (1971). 


1-3 Statement of the Problem 

The present study has evolved from a need to improve on the 
existing methods of the rainfall-runoff relationship analysis, which is 
part of a general problem of hydrograph forecasting. A host of mathe- 
matical models simulating rainfall input and runoff output from a 
watershed are available but there is as yet no general theory to explain 
the course of mechanics of flow from input to output. It is not unusual 
to obtain predictions from different models that differ, to an unacceptable 


degree, both with each other and with the truth. Some of the more obvious 
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causes of this lack of reliability are the following (1., 2., 3. after 
Amorocho, 1967): 

(1) Time variability of watersheds due to man-made changes, and 
due to the natural processes of weathering, erosion, climatic change, etc. 

(2) Uncertainty with respect to the space and time distribution 
of the inputs and outputs of hydrologic systems, and with respect to the 
states and properties of the interior elements of the system in time. 

(3) The inherent nonlinearity of the processes of mass and energy 
transfer that constitute the hydrologic cycle. 

(4) The use of natural data which limits the rainfall-runoff 
analysis to derivation of only empirical and statistical relationships. 
The basic laws or principles involving a certain factor cannot be 


readily determined. 


1-4 Objectives and Scope 


The purpose of the present study has been to extend and improve 
the existing knowledge concerning the relationship between rainfall and 
runoff. The primary objective was to apply state space theory of systems 
analysis to formulate a general mathematical model of watershed runoff. 
In the state space formulation of an n-th order physical system, the 
future performance of the system can be related to the values of n state 
variables which are defined at any given time. The attractiveness of 
the state space method to mathematical simulation in hydrology lies in 
the (1) flexibility and the systematic manner with which the method 
handles any number of inputs and outputs, (2) clear physical interpretation 
of the state variables, and (3) applicability of the method to both 


linear and nonlinear systems. Specific objectives were: 
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(1) To develop and test a state variable model, able to 
simulate surface runoff generated by temporally and spatially distri- 
buted rainfall. 

(2) To compare the state variable model with two existing theories 
of surface runoff: (a) the instantaneous unit hydrograph theory, and 
(b) the kinematic wave theory. 

(3) To investigate the degree of nonlinearity involved in the 
rainfall-runoff relationship, with the aim to assess the feasibility of 
using the theory of linear systems (eg., the unit hydrograph method) in 
the analysis of the rainfall-runoff relationship. 

Because of the weaknesses of the conventional approaches to the 
rainfall-runoff relationship studies (due to the use of natural data), 
the analysis in the present study has been based on the results of 
laboratory catchment experiments. The investigation has been limited 
to rainfall-generated surface sheet runoff from a uniformly sloping 
impervious plane of uniform roughness. These restrictions, which can 
be removed in the future when the mechanics of the surface runoff becomes 
clearer, do not impair the validity of the present approach of 
analysis. On the contrary, the controlled laboratory experiments 
eliminate the uncertainty embedded in natural data, and the non- 
stationarity of natural watersheds. The use of an impervious surface 
eliminates the arbitrariness of the base flow separation. Therefore, 
the concepts and theories concerning the rainfall-runoff relationship 


can be readily examined from the results obtained by this research. 
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CHAPTER II 
LABORATORY CATCHMENT: EXPERIMENTAL EQUIPMENT AND TECHNIQUES 


2-1 Introduction 

. Dimensional analysis is a well established technique in 
engineering fluid mechanics. A useful practical by-product of such 
analysis has been the discovery of how, within practical limits, to 
devise and run models so “rp quantitative predictions can be made about 
prototype behaviour despite total ignorance of the dynamically correct 
formulas for the phenomena at work (Blench, 1969). 

The use of dynamically similar models to solve hydrologic 
problems has been an attractive possibility for many years. In such 
studies the topographic features of a natural (prototype) catchment are 
reproduced to an approximate scale in a laboratory model from which 
predictions of prototype behaviour are obtained by application of 
Similarity criteria. Unfortunately, studies of overland flow conducted 
by Chery (1966), Grace (1967), et. al., show that true hydraulic simili- 
tude between a prototype basin and a scale model is very difficult if 
not impossible to obtain. 

An alternative approach, used in this study and previously 
employed by Amorocho (1961), is to study a laboratory basin, which 
although not representative of any,specific)prototype, still, retains 
some of the qualitative characteristics of natural catchments. By 
analyzing experiments on such basins the aim is to develop a method 


equally well applicable for analysis of natural catchments. 
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2-2 Requirements 


Laboratory study of a small catchment involves a number of 
requirements pertinent not only to the instrumentation necessary to 
control the input and to measure the output, but also to the character- 
istics of the catchment itself. Such requirements are discussed in this 
chapter. Also the designed equipment and the applied experimental 


techniques are described. 


2-2-1 Test Basin 

The size of a laboratory basin should be such that the use of 
simple and reliable instrumentation is possible. The accuracy of the 
outflow measurement increases with the size of the test basin. However, 
with increasing catchment area the difficulty arises in providias 
adequate rainfall input control. 

Depending on the scope of experiments the test basin should pro- 
vide for the possibility to vary the slope, and roughness of the runoff 
surface. In the case of experiments including the subsurface flow, the 
test basin should provide a layer of porous material and allow for 


measurements of the subsurface outflow. 


2-2-2 Input Control and Output Measurement 


The objectives of-using a laboratory catchment demand that any 
spatial variation of input rainfall should be reflected in the runoff 
hydrograph. Thus the ability to produce controlled inputs both in 
time and space is a prerequisite for the detection of such effects. 

The simplest input is a uniformly distributed rainfall in space 


and time. The uniformity of input rainfall is closely related to the 
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method used to simulate rainfall. Basically there are two types of 
rainfall simulators. One type uses sprinklers or nozzles to form rain- 
drops, the other type employs capillary tubes. Sprinklers and nozzles 
produce raindrops of variable sizes comparable to the drop size distri- 
butions observed in natural rainfalls (Hall, 1970). However, it is 
difficult to achieve uniform areal distribution. To maintain the 
uniformity of the inputs and to investigate the effect of raindrop 
impact on sheet flow, it is desirable to limit the experiments to one 
drop size at a time. 

In summary, the input installation should have the following 
features: 

(1) Drop sizes comparable with those occurring in natural 
rainfall. 

(2) Allow variation in rainfall intensity. 

(3) Assure uniform application over the test basin. 

(4) Ability to produce time and space variable inputs. 

(5) Quick response to on-off inflow commands. 

From the existing methods of simulating rainfall (Hall, 1970) 
only the method employing capillary tubes to form raindrops appears to 
satisfy the above requirements. 

The outflow measuring device should be simple, reliable, and 
accurate. It is desirable to measure outflow directly, however, no 
such method, applicable to sheet flow, is available. Instead, the 
cumulative volume or weight of outflowing water is usually measured. 
Relatively simple transducers of the cumulative amount of outflow can 


be constructed by utilizing strain gauges (Byards, 1969; Stein, 1962). 
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2-3 Experimental Equipment 


All devices used in the laboratory catchment experiments, except 
the data acquisition system, were designed and made for the purpose of 
this investigation. A full description of these devices is given in the 


following sections. 


2-3-1 Apparatus 


The general arrangement of the experimental set up is shown 
schematically in Figure 2-1. The photographs of the installation are 


shown in Figures 2-2 and 2-3. 


2-3-2 Rainfall Simulator 

The rainfall simulator, for reasons given earlier, was of the 
type in which drops dre formed by a large number of short length fine- 
bore capillary tubes. Mechanics of drop formation from a capillary tube 
was described by Chow (1965) and need not be discussed here. 

The unit, shown in Figure 2-4, consisted of a closed box, 24 1/4 
inches by 37 1/2 inches, constructed of 3/16 inch thick plexiglass. The 
height of the box increased diagonally from 7/8 inch to 3 1/8 inch. The 
sloping ceiling allowed concentration of air bubbles in one corner of 
the box where a vertical tube served both as an air vent and an open 
manometer. 

The water supply was provided by a 3/4 inch copper pipe, placed 
inside, across the width of the box, at the higher end. Drilled holes 
along the pipe faced the top of the box to minimize the effect of inflow 
velocity. 


The raindrop producing tubes were 1 inch long, made of stainless 
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steel, with 0.012 inch I.D. They were placed at 1 1/4 inch centers in 
the bottom face of the box, (Figures 2-5 and 2-6) totalling to 532 tubes. 

The unit rested on an aluminum rectangular frame provided with 
four levelling screws. The frame was placed on top of an overall 
framework. 

Pressure inside the unit was controlled by a constant head tank, 
shown in Figure 2-7, movable along four vertical rods. Maximum possible 
difference in pressure head, between the water level in the tank and 
the tips of the capillary tubes, was 3 feet. This installation was 
suitable for generation of steady rainfall inputs and a ramp function 
type of inputs. 

It was intended to monitor the pressure inside the unit by a 
pressure transducer, calibrated in terms of rainfall intensity. However, 
the pressure transducer was found to be temperature sensitive, and since 
the daily temperature fluctuations in the laboratory were up to 10°F in 
summer, this idea was abandoned. Instead, rainfall intensity was 
determined from measured steady outflow during experiments. 

The rainfall intensity range of the described simulator was from 
0.75 inches per hour (19.05 mm/hr) to 6 inches per hour (152.4 mm/hr). 
Since the unit consisted of 532 tubes with inside diameter = 0.012 inches 
(0.3048 mm) and covered an area of 6 square feet (5574.1824 cm2), the 


average flow velocity in cm per second through the tubes was 


5574.1824 x i LOTGEST A (2-1) 


ee 03048 
4 


3600 x 532 x 


where i is the rainfall intensity in mm per hour. For the proposed 


range of rainfall intensities, the average flow velocity in the tubes 
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varied from 7.58 cm/sec. to 60.79 cm/sec. The maximum Reynolds number 
(Re = ”)s was thus about 200, well within the range of laminar flow, 
and the flow was therefore assumed to be laminar. This in turn assured 
a linear relationship between applied pressure to the unit and the resulting 
rainfall intensity (Amorocho, 1965). 

The length of the capillary tubes, L, was found to be an important 
factor, greatly influencing the speed with which the individual drops 
are formed. For laminar flow the friction head loss, hg, is expressed 


by the Hagen-Poiseuille equation 


eesZVLV 


where v is the kinematic viscosity of water and g is the acceleration of 
gravity. For maximum velocity of 60.79 cm/sec the head loss per 1 mm 
of a tube's length is hf = 1.98 cm. This represents about 4% of the 
total pressure range used in the experiments (about 50cm), or 4% change 
in the corresponding range of rainfall intensity. It is therefore 
important that the tubes are as much as possible of the same length to 
assure spatially uniform rainfall distribution. 

The average height of the rain simulator above the test basin was 
4.4 feet (134 cm), and the average diameter, d, of a drop was 2.5 m. 
The average drop diameter was estimated from a sequence of photographs 
taken by a movie camera as illustrated in Figure 2-8. By interpolation 
from Wenzel's data (1970), shown in Figure 2-9, giving the relationship 
of drop size, fall height and fall velocity, the impact raindrop 


velocity in the present study was approximately 15 ft/sec. 
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Figure 2-1. 
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Sketch of the laboratory catchment apparatus installation. 
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Figure 2-4. Rainfall simulator with a constant head tank 
and a rain-intercepting tray used to time the duration of 


rainfall applied to the laboratory catchment. 


Figure 2-5. View of the 2 by 3 feet rainfall simulating 
unit showing 532 stainless steel tubes (1.D.=0.012 inch) 


spaced at 1 and 1/4 inch intervals. 
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of raindrops. 


Figure 2-7. Arrangement of 
a movable constant head tank 
supplying water into the 


rainfall simulator. 
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t = 0.0 sec t= rye sec b= 2 Se sec 


t = 16/18 sec 


tu=. 26/18 sec t=) 27/16 /siec 


Figure 2-8. Time sequence of formation of drops from capillary 
stainless steel tubes (I1.D.=0.012 inch) at the rainfall intensity 
of approximately 3 inches per hour. Average maximum drop diameter 
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2=325,liming Device 

Starting and stopping of rainfall application was done manually 
using a rain intercepting tray, shown in Figures 2-2, 2-3, and 2-4. The 
tray, provided with four wheels, was placed between the rain simulator 
and the test basin, on a separate frame with two leading rails. The 
slope of the rails was 0.49, which enabled the released tray to slide 
down fast by the action of gravitational pull alone. To clear or screen 
the test basin took on the average 0.5 second. 

Water intercepted by the tray was discharged elsewhere by a 
plastic tube. This type of rain interception was also suitable for 


simulation of moving storms over the test basin. 


2-3-4 Test Basin 

A rectangular shape of the test basin, with uniformly sloping 
runoff plane, as shown in Figures 2-10 and 2-11, was used exclusively. 
The basin, made of 1/4 inch plexiglass, was supported by an aluminum 
frame allowing variation of the slope of the basin within the range from 
0 to 20 degrees. 

The runoff plane was made of galvanized sheet metal, treated with 
diluted solution of hydrochloric acid. This treatment changed the 
originally non-wetting metal surface into a wetting surface. 

Experiments with non-wetting surfaces showed that water formed on 
the surface as small globules. The globules grew, and when they became 
of sufficient size, they would suddenly leave their place and flow over 
the surface as a small slug of water. As a result, the runoff hydrographs 


were of undesirable shape in the form of cumulative stepwise increments. 
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2-3-5 Discharge Measuring Device 

The choice of the method of measuring the discharge from the test 
basin was primarily dictated by the amount of outflowing water. Since 
the discharge from the test basin was relatively small and occurred as 
sheet flow with shallow depth, the only practical method was to record 
cumulative outflow continuously. The instantaneous rates of discharge 
were then determined by numerical differentiation. 

Discharge from the test basin was received in a painted sheet 
metal trough, supported by two steel beams. The arrangement of the 
trough and beams is shown in Figures 2-12, 2-13, and 2-14. The design 
permitted only one position of the trough with respect to the beams 
(important for calibration), and minimized friction during bending of 
the beams. 

Only one of the beams was active, with two strain gauges attached 
centrally both at the top and the bottom. The strain gauges, type Budd 
C-6-141-1000, gauge factor 2.02 + 0.5%, were cross — connected into a 
full Wheatstone bridge, as shown in Figure 2-15. The input voltage into 
the bridge was 24V. The adopted strain gauge measuring system was based 
on the linear relationship between weight applied to a beam and the 
resulting strain. This method was preferred to others because the use 
of strain gauges permits to construct relatively simple measuring 
devices, with high sensitivity and accuracy. Also calibration of such 
devices does not change appreciably with time, and modern strain gauges 
are temperature compensated. 

The strain-resistance relationship, as well as the criteria 


important for choice of strain gauges are discussed in the Appendix A. 
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2-3-6 Data Acquisition System 

The signal from the Wheatstone bridge was either continuously 
recorded by a Hewlett-Packard 7004A X-Y recorder or punched on a paper 
tape in form of a discrete time series, using Digitec Data Acquisition 
System. 

The X-Y recorder was used only in a few experiments, dealing with 
the effect of raindrop impact on runoff. For the main series of runoff 
experiments the discrete data system was used exclusively. A flow chart 
of the discrete data acquisition system is shown in Figure 2-16. 

The Wheatstone bridge signal was first amplified ten times before 
the data were punched on a paper tape in ASCII code. The fastest 
possible rate of data punching was about one nine-digit reading per 1.6 
sec. This was much slower than the manufacturer's suggested rate of 
reading (personal communication), however, it was observed that with 
increased frequency of reading the system malfunctioned. 

For further processing the data were copied on a nine-track 
magnetic tape, which was subsequently used in IBM-360 digital computer. 
The final output was a printout of cumulative and instantaneous outflows 
as functions of time, and computer plotted hydrographs of instantaneous 
discharge. 

Despite a certain noise level present in the electronic system, 


the load readings on the digital system could be made to within + 3g. 


2-4 Experimental Techniques 


Timing of the rain application to the test basin was done manually 
by releasing or pulling up the intercepting tray. This technique 


inevitably introduced random errors into the measurement. Each experiment 
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was therefore repeated seven times to obtain an estimate of the average 
runoff values with respect to the assumed beginning and end of the rain- 
fall. 

This, in turn, required that the initial conditions, with respect 
to the amount of water present on the runoff surface before each 
repetition, were the same. To start each run with initially dry surface 
would have been too time consuming. Morgali (1970) concluded from his 
analysis of the scatter of the data obtained in laboratory catchment 
experiments, that the ideal condition is for the runoff plane to have 
neither free water on the surface, which would add to the flow, nor 
dryness, which would cause water to be absorbed by the plane. Hence, 
five minutes drying interval was allowed after each run before the next 
run was started. This period was sufficient for the runoff from a 
previous test to cease to zero, within the limits of the recording 
device. 

Calibration of the outflow measuring system was done for weight 
increments of 500 grams by adding 500 cc of water into the collecting 
trough. The relationship between applied weight and the Wheatstone 
bridge voltage output, shown in Figure 2-17, was found to be linear 
within the working range from 0. to 11,000. grams. The calibration was 
checked before and after each experimental period, and was found to be 
stable. 

The noise in the recorded signal was partly due to the fluctuations 
of the reference voltage applied to the Wheatstone bridge. The fluctu- 
ations were observed to be within the range 24.00 + 0.02 volts. Another 
source of noise was a slight vibration of the beams caused by the impact 


of inflowing water into the trough. 
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The experimental runs were conducted in the following manner: 

(1) Approximately 20 minutes warm up period preceded the beginning 
of an experimental period. During this period the electronic parts of 
the system were allowed to stabilize. 

(2) The constant head tank was secured in position for a desired 
rainfall intensity, while the sliding tray was in position intercepting 
the rain. 

(3) After the warm up period the sliding tray was released and 
the test basin was wetted completely. 

(4) When a steady outflow was reached the rain was interrupted. 
After a five minutes drying period the first experimental run started. 
The initial voltmeter reading was recorded and an electric stop watch 
and the paper tape punch were started. Simultaneously the sliding tray 
was released to begin the rain application. 

(5) After a specified duration of rain, the sliding tray was 
pulled up and the receding runoff continued to be recorded for another 


150 seconds. 
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Figure 2- 


runoff 


10. Laboratory catchment (size 2 by 3 feet) and 


collecting trough resting on two steel beams. 


Figure 


2-11. Sheet flow over the laboratory catchment 


runoff plane under simulated rainfall. 


30 


or) 


ny oe “v4 = ate te 


wre (laod € vd S eade) dae varosbs Yrm ako LOE. 
1 a . . / dari 


Lempad Iegte. owt ate irks tate | iguord aiktael too, Hames | 


DM {i 
xh } hae 
_ ; rn 
4 
; <wih 
. hy 
‘4 4 i 
f Be em 
A yy 
i iy yi a 
' 
eat 


a aie j om 


en i 


31 


°ueeq pecoT ey} FO YOISYS °7ZI-Z sansty 


xi xv 
Wv4a ONILYOdaNs - 


8311048 
aI ONIWV3S8 11V98 


pie 


1yOddns \ 
39a3 34INX 


WvV3d JONVS Nivals NO1d31 


g-q NOILDas 


GOxGxG 
ANVL ONILDIIVIOD 


W- NOILD3S 


oN y 


N | a 310H 
4 | Iv DINOD 
PAS | 8 ran 


fe 


; i: i f r) ) ; 
{ ’ F Lea 
Nes 


ty Lame 


; i oi el GO he i Ny iy vag Wt ha ah 
| ‘ifn fy iy it Sate: LL a Mt : jy ; eal 
a y y beak ely i MON) lB hl ind ni 
( mare ‘ Tae A i i 4 , ne : F 
,) fi an 7 ‘ i: i \ , a wa 4 , 
; } FY co ATMA ic pel ; a Aa hh ae ; ‘ 
’ uy " a ; f ak | By [ \ i Leer 
| i ; Te Ve ane Ay 


es~ 


evercp St sys jyosgq p 


= 


ZiBVIM GYNCE2 


a - 


ties 
an 
Rp 
: aG 
ihe Be 
a) fi 
E j 
ve 
4 
) 
‘ 
j : \ . 
i ‘ 
: 4 


em — Geert he Nie ue 
: hie Wey Caen Gi a, ee ey iri Pee 
TOW ane RNS Weae is aut, da 


‘ arn) 


Figure 2-13. 


Installation of the strain gauge transducer. 


Figure 2-14. 


Arrangement of the collecting trough and the 


load beam. 
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Figure 2-15. Load beam - arrangement of strain gauges into a full 


Wheatstone bridge, and the bridge circuit diagram. 
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Figure 2-16, Flowchart of the data acquisition system. 
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CHAPTER III 
NONLINEARITY PROBLEM IN HYDROLOGIC SIMULATION 


3-1 Introduction 

The nonlinear nature of watersheds behaviour has been observed 
and reported in recent years by many authors (Minshall, 1960; Amorocho, 
1961, 1971; Diskin, 1964; Kulandaiswamy, 1964; Proceedings of the 
International Hydrology Symposium, Fort Collins, Colorado, 19673; Chiu, 
1970; First U.S.-Japan Seminar in Hydrology, 1971). Because there is 
no general theory of analysis applicable to nonlinear systems, 
"identification" of a nonlinear system proves to be a formidable task in 
most cases. Because of the theoretical and practical difficulties 
associated with the treatment of nonlinear systems most hydrologic 
analysis have been directed towards the formulation of procedures for 
the linear approximation of these systems (Amorocho, 1967). How close 
such approximations may be is a problem which has not been given much 
attention in hydrology. 

In the following sections the risk associated with the use of 
linear models approximating nonlinear systems is discussed in 


qualitative terms. 


3-2 Linearized Hydrologic System 


Three basic equations are available for the deterministic 
simulation of the relationship between rainfall and runoff. These are: 
(1) equation of continuity, (2) momentum equation, and (3) energy 


equation. When the operation of the rainfall-runoff process of the 
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hydrologic cycle is approximated by a linear model, it is assumed that 
the momentum equation is linear. As a very simple example the effect 
of linearization in the case of a single nonlinear reservoir, as shown 
in Figure 3-1, can be used as an illustration. The single reservoir 
system is similar in some or all respects to other hydrologic-hydraulic 
systems for which the solutions involve relation between storage and 
discharge, such as small dams, ponds, overland flow, and individual 
reservoirs in the reservoir cascade approach to watershed hydrograph 
models. 

The equation of continuity for the single reservoir is 

ds 


Beat) a=) at) (3-1) 


where S is the amount of stored water in the reservoir at time t, i and 
q are time varying inflow and outflow rates respectively. 

The dynamic equation can be written as 

q(t) = K’S (3-2) 
where K = cB 2g/A is a constant; c is the discharge coefficient, B is the 
area of discharge orifice, g is the acceleration of gravity, and A is 
the area of the reservoir. 

Substituting Equation (3-2) into Equation (3-1) gives 


eae oe ‘ 


Response of the nonlinear reservoir to an arbitrary input i(t) can be 
found by solving Equation (3-3). It is interesting to compare responses 
of the nonlinear reservoir and its linear model to a simple input such 
as a step function input. 

Linear model of the reservoir requires that the q-S relationship 


be linearized.Mathematically, linearization of the q-S relation may be done 
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Figure 3-1. 


Definition sketch of a reservoir. 
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by means of a Taylor series expansion of the function q(s), given by 


Equation (3-2), around the steady state value Sy; thus 


q(S) = q(Sx) + q' (Sx) (S-Sy) + WG Vita: (3-4) 
where q'(Sx%) is the first derivative of q evaluated at Sx, q''(Sx) is the 
second derivative, etc. Truncating the series after the last linear 
term, the result is 

d= qx + qx(S - Sx) (3-5) 
where q = q(S), and qx = q(Sx). The first derivative of q, is 


il 
5 


qx = 'KS,°, which after substituting into Equation (3-5) gives 


Bite ol (3-6) 
where the parameter tT -2 s2 is called the time constant of the system 
and has the units of time. (Physically this parameter represents the 

time elapsed since the beginning of a step input applied to a first order 
system operating at a steady state, until the system's response reaches 
63.2 percent of its ultimate value. When the time elapsed is 21, 31, 

and 4t, the percent response is 86.5, 95, and 98, respectively (Coughanowr, 
1965).) Substituting Equation (3-6) into Equation (3-3) yields 


dS ] Dim De ce 
dt fe ane yee oo 


Equation (3-7) is a linearized version of the nonlinear relationship 
given by Equation (3-3) and represents a tangent line passing through 
the point of operation O(q,,Sx), as illustrated in Figure 3-2. Solution 
of Equation (3-7) for a step input of magnitude I (Coughanowr, 1965) is 


S(t) = S, + tI(1 - e7t/T) (3-8) 
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Figure 3-2. Graphical interpretation of linearization of a function 
around a steady state operating point, O(q*,S*), by means of a 


Taylor series expansion. 
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3-2-4 Comparison of Measured and Computed Responses 


Predicted responses by the linear model to step function inputs, 
were compared to the true nonlinear responses determined by measuring 
the changes of storage in a small plexiglass tank whose nonlinear 
function of outflow was given by Equation (3-2) with K = 0.74. The tank 
was placed on two beams, one of which was fitted with strain gauges, 
thus acting as a weight transducer (Chapter II). The signal from the 
transducer was recorded in time intervals of 4 seconds. Outflow from 
the tank occurred through a circular orifice centered in the bottom of 
the tank, and the inflow was regulated by a valve. 

The experimental measurements consisted of two step inputs 
applied to the tank initially operating at a steady state (same for 
both inputs) given by: Sx, = 1955 cm3, ix = qx = 33 cm3/sec. The 
corresponding value of the time constant of the linear model was 
t = 119.5 sec. In the first experiment the inflow was suddenly 
increased to 35.6 cm3/sec, which represents 7.9 percent of the initial 
value of inflow ix. In the second experiment the inflow was increased 
to 137.4 cm3/sec, or by 13.4 percent of the initial value. Measured 
responses to these two step inputs, and the predicted responses by the 
linear model are shown in Figure 3-3. At time t = 4t the predicted 
responses are 13 percent and 19.4 percent respectively less than the 
true values, in terms of the measured increase in storage from the time 
Cs Ustoot = 47, 

The results show that even when the inputs varied within a 
relatively small range (less than 14 percent) around the operating point, 


linearization of the system's equation resulted in the loss of ability 
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to predict responses of the system quantitatively. Similar situations 
arise when dealing with the rainfall-runoff relationship. It is shown 
in Chapter IV that the instantaneous unit hydrograph derived for the 
laboratory catchment is a linear model based on the rainfall-runoff 
relationship linearized for a particular rainfall intensity input. When 
the rainfall intensity changes the particular instantaneous unit hydro- 


graph becomes a poor model of the rainfall-runoff process. 
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CHAPTER IV 


EXISTING MODELS OF SURFACE RUNOFF 


4-1 Introduction 

In the study of surface runoff from watersheds numerous attempts 
have been made to formulate the relationship that is required for the 
simulation of overland flow. These efforts have resulted in the 
development of a multitude of mathematical models describing the runoff 
phenomenon(International Hydrology Symposium, Fort Collins, Colorado, 1967). 

The models can be classified into two categories. The first of 
these categories covers the runoff models which assume a linear relation- 
ship between rainfall and runoff. Most of the linear models are based 
on the unit hydrograph theory. In the second category fall the models 
which are nonlinear. 

In the present chapter two of the more popular concepts of surface 
runoff are reviewed and performance of mathematical models based on 
these concepts was tested experimentally. These two concepts are the 
instantaneous unit hydrograph and kinematic wave concepts. 

The IUH model was chosen to represent all linear models of the 
surface runoff in this study, primarily because the unit hydrograph 
theory has become very familiar to most hydrologists. If the investiga- 
tion reveals that the nonlinearity of the surface runoff is a significant 
reason for the IUH failure, then it can be concluded that any linear 


model is susceptible to the same degree of failure. The kinematic wave 
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model is a nonlinear model based on the kinematic wave theory. The 
attractiveness of the model lies in the relatively easy solution of a 
simplified hydrodynamic equation governing the runoff. 

The unit hydrograph and the kinematic wave have been considered 
theoretically for many years. Both concepts have been used quite 
extensively in attempts to solve practical problems. However, there 
seems to be a lack of experimental investigation of their validity and 
properties. Review of the requirements phdeccndrtaadk assumed in the 
unit hydrograph and kinematic wave theory indicates that these concepts 


can be best tested by means of a laboratory catchment. 


4-2 The Instantaneous Unit Hydrograph 


The instantaneous unit hydrograph (IUH), designated in the 
following by u(t), is a unit hydrograph resulting from an effective 
precipitation of infinitesimally small duration. In other words, for an 
IUH the effective precipitation is applied to the drainage basin in zero 
time, and the input is in the form of a unit impulse. This rather 
fictitious hydrograph has a major advantage in comparison with a unit 
hydrograph. It is independent of the duration of effective precipitation, 
thereby eliminating one variable in hydrograph analysis. Furthermore, 
the IUH is better suited for the needs of theoretical investigations of 
the rainfall and runoff relationship because it allows the use of the 


linear theory of systems analysis. 


4-2-1 The Convolution Integral 
It is known from the systems analysis (Zadeh, 1963; Bendat, 1966) 


that the dynamic characteristics of a constant parameter linear system 
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can be described by a unit impulse response, h(t), defined as the output 
of the system at any time t to a unit impulse input applied at time t-t. 
If a catchment is considered to be a linear system with the surface run- 
off as an output, then the IUH, or u(t), is the system's unit impulse 
response and the effective rainfall of the IUH corresponds to the unit 
impulse input. 

Usefulness of the unit impulse response for description of the 
system is due to the fact that for any arbitrary input i(t), the system's 
output Q(t) at any time t is given by the convolution integral 

c 
Q(t) = | i(t) u(t-t) dt (4-1) 
Oo 
According to Equation (4-1) the value of the output Q(t) at time t is 
given as a weighted linear sum over the entire history of the input i(t). 
Once the unit impulse response, sometimes called the weighting function 
of the system, is known the system is "identified", and the system's 
output can be predicted for any input by Equation (4-1). The prediction 
property of the convolution integral led to its use in hydrograph fore- 
casting. The assumption of linearity of the surface runoff process was 
essential. In this approach the problem of runoff simulation was then 


reduced to finding the unit impulse response of a catchment. 


4-2-2 Conceptual Models of IUH 


There have been many conceptual models proposed to delineate an 
IUH. They were composed of simulated components such as linear reservoirs 
(Nash, 1960), linear channels (Dooge, 1959) or time area diagrams (Ross, 


1921). The theoretical development of the IUH models composed of linear 
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reservoirs has been derived in the following sections using the transfer 
function method of linear systems theory. The method, as well as the 
conclusions about the reservoir type model, generally apply to any linear 
model. The reservoir type model was chosen merely because it can be 
compared with the nonlinear state variable runoff model derived in 


Chapter V, and which is also composed of reservoirs. 


4-2-2-1 Single Reservoir 


The simplest conceptual model of a catchment is a fictitious 
reservoir in which the outflow q is directly proportional to the storage 
S(Amorocho, 1961), or 

S = kq (4-2) 
where k is constant for the reservoir and has units of time. The inflow 
i into the reservoir and the outflow q are related to the rate of change 


of storage by the continuity equation 


dS -i-q (4-3) 


The unit impulse response, u(t), of the reservoir (Appendix B) is 


ait) =< ent/k (4-4) 


Substituting Equation (4-4) into the convolution integral, Equation (4-1), 
and assuming a constant input i, gives the outflow function as 

q(t) = i(1 - et/k) (4-5) 

If the inflow terminates at time tg since the outflow began, then 

q(t) = qoe7*/k (4-6) 
where tT = t - tg, and qo is the outflow at time ty. 

The disadvantage of a one stage first-order system, such as the 


single reservoir, in runoff simulation is that the rate of change of the 
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Figure 4-1. Characteristic responses of non-interacting 


first-order linear systems to unit step function input. 
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response (slope of response curve) is maximal at t = 0. In order to 
compensate for this, conceptual models composed of a series of reservoirs 
have been suggested. Step response of a system consisting of two or 

more first-order systems is S-shaped and the response changes very slowly 
just after introduction of the step input. This sluggishness or delay 

is sometimes called transfer lag (Coughanowr, 1965), and is always 
present when two or more first-order systems are connected in series. 
Figure 4-1 illustrates the step input response curves for several 


systems composed of one or more linear reservoirs in series. 


4-2-2-2 Series of Reservoirs 

Nash (1960) made a significant contribution to the unit hydrograph 
theory by presenting a conceptual model composed of a series of n 
identical linear non-interacting reservoirs, as shown in Figure 4-2. 
All reservoirs have the same storage characteristic, S = kq. Response 
of such a system to a unit impulse input (Appendix B) is 


=I 
Qn(t) = nae or cali (4-7) 


where Once) is a deviation variable of the outflow from the n-th reservoir 
at time t. Nash substituted the factorial in Equation (4-7) by a Gamma 
function to account for continuous values of n. Thus the Nash's equation 


of the IUH is 


Bae wt oc silat t/k 
COE every a) we (4-8) 


Values of k and n can be evaluated by the method of moments from the 


knowledge of past hydrographs. 


There are several basic objections to the Nash model and the unit 


hydrograph concept generally. Firstly, there is the problem of assessing 
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Figure 4-2, Illustration of the physical analog of the Nash 
IUH model as the response, qn, of a series of n non-interacting 


linear reservoirs to a unit impuls input i(t). 
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the correct amount of the effective rainfall as an input into the system. 
Secondly, the IUH method is a black box type of analysis, and as such it 
disregards any information about the mechanics of water motion. A good 
fit is the sole objective of the analysis and the only physical argument 
may be that the runoff process is a heavily damped linear system. 

The Nash model does not provide for spatial distribution of rain- 
fall since all the input into the system enters only through the first 
reservoir. Another criticism that can be made is that the model allows 
for reservoir storage effects but not for channel translation effects 
present in any catchment (O'Donnell, 1966). 

The spatial distribution of the input rainfall and the channel 
translation effects were subjects of later proposed IUH models (Dooge, 
1959; Singh, 1964; Diskin, 1964). However, the single main drawback of 
all IUH models remains the assumption of linearity as can be seen from 


the experimental results presented in the following sections. 


4-2-3 Direct Determination of IUH 

Deriving a catchment IUH by invoking the aid of some catchment 
model is bound to lead to an approximation of the real IUH, the goodness 
of the approximation being dictated by the goodness of the model 
simulating reality. One would like to use methods of deriving an IUH 
that by-pass the need for a model, i.e. operate directly on the rainfall 
excess and direct runoff data to yield the IUH. Such methods have, in 
fact, been used (O'Donnell, 1966) and require only that the system be 
linear and time-invariant. 

The method of deriving the IUH used in the present study was first 


proposed by Chow (1962), It is based on the fact that an S-hydrograph is 
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an integral curve to the IUH. The IUH ordinate at time t is thus simply 


equal to the slope of the S-hydrograph at time t. 


4-2-4 Comparison of Experimental and Theoretical Results 


The experimental apparatus, described in Chapter 2, was well 
suited for experimental testing of the unit hydrograph theory because 
the two most basic assumptions of the theory were satisfied: (1) the 
input rainfall was uniformly distributed over the experimental basin, 
and (2) the measured hydrographs resulted from direct runoff only. 

The aim of the experiments was not to test particularly any of 
the IUH conceptual models. Rather, the objective was more general: to 
depict clearly nonlinearity of the runoff process and thus to prove 
invalidity of the unit hydrograph method generally. This amounted 
basically to: 

(1) Derivation of an IUH from a recorded experimental S-hydro- 
graph caused by a rainfall of specific intensity. 

(2) Recording a number of experimental hydrographs resulting from 
rainfalls of different intensities and durations. 

(3) Simulation of the experimental hydrographs by means of the 
convolution integral, using the derived IUH as the unit impulse response 
of the laboratory catchment. 


(4) Comparison of the experimental and the simulated hydrographs. 


4-2-4-] IUH of the Laboratory Catchment 


The method of direct determination of the IUH described in Section 
2-3 was used to derive the IUH of the laboratory catchment. The method 
requires that a continuous record of cumulative outflow from the catchment, 


i.e. S-hydrograph, be available. An S-hydrograph was easily obtained 
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experimentally by letting the duration of rainfall be longer than the 
time to equilibrium. A smooth curve was then fitted to the experimental 
data and the IUH was obtained by numerical differentiation of the curve. 
Two recorded S-hydrographs, S; and S9, resulting from 3.84 cm/hr 
and 11.53cm/hr rainfall intensities, and the corresponding derived IUH's, 
IUH-1 and IUH-2, are shown in Figure 4-3. The two IUH's were adjusted 
to the common "unit" rainfall intensity of 10 cm/hr. A linear rainfall - 
runoff process would have resulted in derivation of exactly identical 
IUH's from both S-hydrographs. The actually derived IUH's show 
remarkable differences in peak values, time to peak, and the total 
duration. Hence, the sheet runoff caused by rainfall must be classified 


as a highly nonlinear process. 


4—-2—4-2 Comparison of Simulated and Observed Hydrographs 


Hydrographs simulating recorded runoffs from the laboratory 
catchment, generated by rainfalls of various intensities, were calculated 
by convolution of the derived IUH's with the input rainfall intensities. 
Since the experimental data were discrete, summation form of the 
convolution integral (Chow, 1964) was used, expressed as 


n 
Oe) = 2 ult - G = 1)dt}ijat (4-9) 
jai 


where t represents any time after the beginning of effective rainfall. 
The input rainfall may be considered to consist of n blocks of different 
intensities I; and of the same duration t; the subscript j denotes the 
number representing a block. A computer program performing the 


convolution given by Equation (4-9) is listed in Appendix D. 


The convolution results provided an interesting insight into the 
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Figure 4-3. Recorded S-hydrographs, Sj and S9, of runoff from 
the laboratory catchment, and the corresponding derived instan- 


taneous unit hydrographs IUH-1 and IUH-2. 
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IUH simulation method. Typical results are shown in Figures 4-4 and 
4-5. Three experimental hydrographs resulting from the rainfall input 
of the same intensity as for the Sj-hydrograph were simulated. The 
rainfall input durations for the three hydrographs were 20, 25, and 50 
seconds. 

Firstly, the hydrographs were simulated using the IUH-1 as a unit 
impulse response in the convolution equation. To repeat the earlier 
comment, this was actually a simulation of hydrographs caused by the 
same rainfall intensity as was the Sj-hydrograph from which the IUH-1 
was derived. Figure 4-4 shows that the rising parts of the hydrographs 
were reproduced exactly up to the time equal to rainfall durations. 
However, the recession curves were extremely poorly simulated and the 
times to peak and the peak values were grossly overestimated. 

The overshoot of the simulated hydrographs is due to the shape of 
the IUH-1. The IUH-1 peaks approximately at t = 19 sec, as can be seen 
in Figure 4-3. Thus the past inputs are weighted more heavily than an 
immediate input. 

Secondly, the same three hydrographs were simulated using the 
IUH-2 as a unit impulse response in the convolution. Figure 4-5 shows 
the results which strikingly demonstrate nonlinearity of the runoff pro- 
cess. The simulated hydrographs in this case do not even partly resemble 


the experimental curves. 


4~2-5 Nonlinearity of Runoff Distribution 


Comparison between measured hydrographs and those predicted by the 
IUH method indicates clearly that the runoff phenomenon cannot be analyzed 


successfully by linear methods. To illustrate further the nonlinearity 
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of the rainfall-runoff relationship the following experiment was per- 
formed. The runoff plane of the laboratory catchment was fixed in the 
position such that the bed slope was S, = 0.2079. Hydrographs of surface 
runoff from the catchment were recorded for the following three rainfall 
intensities! 3.904 é¢m/hry..7./7 -cin/hit, atid 11.53 .cm/hr. The Weineali 
intensities were approximately in ratios 1:2:3. Duration of rainfall 
for all events was 50 seconds. Average observed hydrographs for each 
rainfall intensity are shown in Figure 4-6. In the same figure are 
plotted two hydrographs, corresponding to rainfall intensities of 

7.77 cm/hr and 11.53 cm/hr respectively, derived by superposition of the 
ordinates of the recorded hydrograph generated by rainfall of 3.84 cm/hr 
intensity. These two derived hydrographs represent in effect a linear 
relationship between rainfall input and runoff output. The actual 
rainfall-runoff relationship, however, is highly nonlinear as can be 
seen from comparison of the recorded hydrographs with those derived by 


superposition. 


4-2-5-1 Time Distribution of Runoff 

In design flood estimation characteristic times are required for 
the determination of hydrograph parameters and critical durations of 
flood-producing rainfall. The rise time, time of concentration, time to 
equilibrium and volume/peak ratio have all been suggested as appropriate 
characteristic times (Bell, 1969). These are influenced to a great degree 
by the physical configurations of the catchments concerned. Many of the 
suggested formulas in fact define the characteristic response times only 
in terms of catchment parameters, such as the slope and length (Snyder, 


1938; Kirpich, 1940; Taylor, 1952). The unit hydrograph theory in effect 
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also assumes a constant response time. However, because the rainfall- 
runoff relationship is nonlinear, the time distribution of runoff is a 
function of the input amplitude, i.e. the rainfall intensity, besides the 
catchment parameters. 

To support this reasoning by experimental evidence the 
relationship of the apparent time to equilibrium, rainfall intensity and 
the runoff surface slope was investigated experimentally. Five different 
slopes of the runoff surface were tested. For each slope the apparent 
times to equilibrium were measured for between 23 and 26 different rain- 
fall intensities within the range from 1.7 cm/hr to 14 cm/hr. The 
results are given in Table 4-1. The apparent times to equilibrium were 
determined from cumulative curves of runoff plotted by an X-Y recorder. 
Typical cumulative curves of runoff are given in Figure 4-7. A tangent 
was drawn to each cumulative curve at the point of equilibrium. The 
intercept of this tangent on the abscissa was taken as the apparent time 
to equilibrium. 

Analysis of the experiments show that the relationship between the 
apparent time to equilibrium, te, and the rainfall intensity, i, can be 
expressed in the form of an equation of a hyperbola, i.e., in general 
by the equation 

tare mee = (4-10) 
where the coefficient a is a function of the slope S,. The effect of the 
length of the runoff plane was not investigated. The experimental data 
for each slope were fitted by least square method with best fit lines 
and plotted on logarithmic paper as straight lines, as shown in 


Figure 4-8. Equations of these lines are also given in the same figure. 
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APPARENT TIME TO EQUILIBRIUM, ta, AS FUNCTION OF RAINFALL INTENSITY, 


i, AND BED SLOPE FOR A RUNOFF PLANE OF LENGTH 91.4 cm 
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The significant result of these experiments is that the time 
distribution of overland flow is a function of rainfall intensity, in 
addition to slope, surface roughness, length, etc. Hence, any relation- 
ship for catchment characteristic times, such as "time of concentration", 
"basin lag", "time to equilibrium", and "time-area concentration curve 
(Nash, 1958)'', which all have been used as measures of the speed of a 
catchment response to rainfall, in fact depends on the rainfall intensity. 
Many of the published relationships (Snyder, 1938; Taylor, 1952; Kerby, 


1959; reviewed by Bell, 1969) failed to recognize this. 
4-3 The Kinematic Wave Model 


4-3-1 Kinematic Waves 

Certain properties of a distinct type of wave motion, normally 
referred to as kinematic waves, were known already at the turn of this 
century. The velocity, c, measured relative to the bank of such waves 
was observed by Seddon (1900) to be given by the slope of the discharge - 
area curve, i.e., 


dQ _ iQ 


fenddAsonBady hae}. 


where Q is the discharge, A is the cross-sectional area of flow, B is 
the channel width, and y is the depth of flow. The validity of Equation 
(4-11), known as Seddon's Law (Linsley, 1958), was demonstrated by Seddon 
on the Mississippi River. 

However, the first proof of the existence of the kinematic wave 
was given by Lighthill and Whitham (1955), who derived its properties 
mathematically. The authors provided the theory for this type of wave 


motion, which arises in any one-dimensional flow problem when there is an 
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approximate functional relation at each point between the flow Q 
(quantity passing a given point in unit time) and concentration k 
(quantity per unit distance). The wave property then follows directly 
from the equation of continuity satisfied by Q and k. In view of this, 
the authors suggested calling these waves "kinematic" as distinct from 
the classical wave motions, which depend also on Newton's second law of 
motion and are therefore called dynamic. 

The kinematic wave theory is based on the assumption that Q is a 
function of y only. Then the total derivative, dQ/dy, in Equation (4-11) 
has a unique meaning derived from the uniform flow condition, i.e., that 
the friction slope Sf and the bed slope S, are equal (Henderson, 1969). 
The assumption that the discharge is a function of depth only is most 
likely satisfied when So is large. Thus, the kinematic wave theory 
appears to be particularly suited for the surface runoff problem, because 
of the steep slopes in the upper catchment regions where overland flow 


exists. 


4-3-2 Kinematic Wave Theory of Overland Flow 


Henderson and Wooding (1964) applied the kinematic wave theory to 
an elementary case of surface runoff illustrated in Figure 4-9(a). The 
runoff occurs from an impermeable uniformly sloping plane, and is caused 
by a uniformly distributed rainfall of constant intensity i. Following 
the derivation given by Henderson and Wooding, the equation of continuity 


for unsteady spatially varied flow caused by rainfall is 


ee ee (4-12) 
ox ot 


where q is the discharge per unit width, y is the normal depth of flow, 
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x is the distance in the main direction of flow, t is the time, and i is 


the rainfall intensity. Equation (4-12) can be written as 


dq oy 4 oy = 7 ues 
dy ox ot (4-13) 
or 
oy eee es 
“ox 5 dt z, (4-14) 
where 
c= Gy for constant x (4-15) 


Equation (4-15) is the Seddon Law and expresses the assumption that the 
discharge at a fixed distance x depends only on the depth at x. 


Since y = fn(x,t) the total derivative of y is 
dy = dx oy 4 dy (4-16) 
Comparing Equation (4-14) and Equation (4-16) one can see that if 


oe aoe (4-17) 


he ee, (4-18) 


Physical interpretation of this result is that, to an observer moving 
with velocity oa = c, the relation given by Equation (4-18) appears to 
be true, as well as the following relations: 

(1) y=it (4-19) 


there being no constant of integration if the surface is initially dry, 
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dq dy _ dg _ . 
(2) a5 deme ci (4-20) 
and 
dqd-d0 dq 2x, 
dt-dx'. dx ©e (4-21) 
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hence 
qi = six + -constant (4-22) 
As a result of the assumption that S¢ = S, the discharge q is 
related to the depth y by the uniform flow equation of the general form 
aeeoye 


S 
Boe m = 3 for laminar flow, and ao = C1Soe, m = 3/2 for turbulent 


where a = 
flow (Chezy equation). v is the kinematic viscosity of water, and g is 
the acceleration of gravity. (For details see Chapter 5, Section 3-2-2.) 
When a steady state is reached, the flow profile is given by 
q = ix (4-23) 
Henderson (1964, 1969) discussed further solution in terms of a 
method of characteristics, i.e., tracing of the wave motion on the x-t 
plane. The resulting curve, called a characteristic, can be visualized 
as a path traced on the x-t plane by an observer moving with velocity c. 
Tracing the movements of many such observers, all starting to move at 
time t = 0 but from different starting positions along the runoff plane, 
will result in a set of characteristics. There is only one set of 
characteristics on the x-t plane, since the kinematic wave motion proceeds 
only in one direction, i.e. downstream. 
Corresponding to each characteristic there is a curve in the 
physical y-x plane. Comparison of Equation (4-22) and Equation (4-23) 
shows that all these curves are identical with the steady state profile 
but displaced from it by an amount depending on the starting position of the 
observer at t = 0. By joining up points on these curves having the same 


value of t, one obtains an instantaneous profile of the water surface 


appropriate to that value of t, as illustrated in Figure 4-9(b). 
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Figure 4-9. Kinematic wave solution of surface runoff by the method 
of characteristics: (a) definition sketch, (b) development of flow 
profiles: AM|N] at time t ;, AM2N2 at time tz, steady state profile 
AM, M5M3 at time t,. Broken lines represent y-x curves corresponding 
to the paths traced by imaginary observers who started to move along 
the plane AB at time t = 0 with velocity c, from their starting 


positions 1, 2, 3, and 4. 
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Since Equation (4-19) appears true to all observers it follows 
that the successive profiles MN in Figure 4-9(b) are parallel to the 
surface AB in the region to the right of the steady state profile. The 
plateau MN rises with constant velocity i until a steady state is 


reached after a time 


mest. L 1/m 4 
io Pam prawn Neer ie (4-24) 


Discharge q at the lower end, B, of the runoff plane is given by 

aimem (t<ts) | 

| Beet re { i (425) 

Li (t>ts) 

The relation for discharge given by Equation (4-25) is valid for 
the duration of rainfall, tj, equal to infinity. For finite durations 
the following two cases of discharge subsidence, after cessation of 
rainfall can be distinguished (derivation of formulas can be found in 
Henderson [1964, 1969]): 

(1) tg2ts. The initial profile is a steady state profile. The 
depth y at B at a given time, tg, after the cessation of rainfall can be 
computed from 


L = aym-1 CF + mtq) (4-26) 


and the discharge is given by q = ay™, 

(2) to<tg. The greatest depth at B reached during the rising 
period is yo = ity. This depth remains constant during the recession for 
the time ta<tp, where tp is given by 

m 
ayo 


L=- ¢ 1 
Loh os tel - attos - to) (4-27) 


L/ayst, From Equation (4-24) it follows that tog = tg when 
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Figure 4-10. Outflow q - t curves. 
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to = tg. When ta>tp, Equation (4-26) holds and the outflow is then 
identical with that given Ey case (1) with the same value of i. 

In summary, three types of hydrographs can be distinguished when 
the surface runoff is simulated by the kinematic wave model. The 
criterion governing the occurrence of a particular type is given by the 
relation of the rainfall duration, to, to the time to equilibrium, tg, 
i.e. whether tp is greater, equal or less than tg. Typical hydrographs 
for each case are shown schematically in Figure 4-10. A characteristic 
property of the kinematic wave model is the sharp crest of the simulated 
hydrograph. This property follows from Equation (4-19), according to 
which the depth of flow increases with constant velocity until a steady 


State is reached or until the rain input ceases. 


4-3-3 Comparison of Simulated and Observed Hydrographs 


In Section 4-3-2 the concept of kinematic wave was used to 
develop the relationships given by Equations (4-24) to (4-27), specifying 
the temporal and spatial distribution of discharge from a uniformly 
sloping plane under rainfall. These equations, which together form the 
kinematic wave model of pera ciao. were used for development of a 
computer program for performing digital simulation of experimental hydro- 
graphs. The computer program is listed in Appendix D. The simulation 
included six events of runoff, from the laboratory catchment with the 
bed slope So = 0.20791, generated by rainfalls of intensities of 11.53 cm/hr 
and 7.77 cm/hr (runs number 1 to 6 given in Table 5-2, Chapter 5). 

The case when the rainfall duration is longer than the time to 


equilibrium, i.e. to>ts, is shown in Figures 4-11 and 4-14. The simulated 
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hydrographs fit the experimental data quite well except near the hydro- 
graph crest. Larger discrepancies in peak discharge values are apparent 
for the case when tp = tg, shown in Figures 4-12 and 4-15 and the case 
when to<tg, shown in Figures 4-13 and 4-16. In these two cases, i.e. 
when to<ts, the predicted peak discharges are considerably higher than 
the observed ones. There are two possible reasons for these discrepancies, 
besides the errors in measurement. Firstly, it is a property of the 
kinematic wave model, given by Equation (4-19), that a steady rate of 
increase of flow depth at a given cross-section occurs until a steady 
state of discharge at the same cross-section is reached. This property 
results in sharp crested hydrographs, contrary to the assymptotical 
approach to a steady flow in natural hydrographs. Consequently, the 
simulated values of discharge near the peak are overestimated. Secondly, 
the raindrop impact effect was neglected in the kinematic wave model. 

As a result, the simulated hydrographs are rising faster and the peak 
discharges are again overestimated. Complete discussion of the raindrop 
impact effect is given in Chapter 5, Section 4. 

The duration ts, sometimes called the "time to equilibrium", appears, 
in one form or another, in all standard hydrological theories or methods. 
In many cases it is assumed to be independent of i, as it was discussed 
in Section 4-2-5-1 of this Chapter. It was found experimentally in the 
Same section that the apparent time to equilibrium, te, is in fact a 
function of rainfall intensity i, and has the form given by Equation (4-10). 
The empirical equations, given in Figure 4-8, specifying te for different 
values of So as a function of i are compared in Table 4-2 to the time to 


equilibrium, ts, predicted by the kinematic wave model. The time to 
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equilibrium, ts, is given by Equation (4-24) which may be written for 
the case of laminar flow from the laboratory catchment (L = 3 feet, 


m= 3) as 


1/3 


1/m 
3Lv 1 Li 
a (ne i yey) im ea ae hae es 
ts “at ? %) 9 © Aaa (4-28) 


The predicted values of tg compare closely with the empirical values of 
te. Because te is of a shorter duration than ts by definition the con- 
stants in the empirical equations are generally smaller than the constants 
in the corresponding equations for tg. The relationship specifying the 
time to equilibrium, ts, follows from the assumption of the kinematic 
wave theory that the discharge is a function of the flow depth only. The 
results in Table 4-2 support the validity of this assumption or at least 


give the evidence of a good approximation for the given conditions. 


TABLE 4-2 
COMPARISON OF THE EMPIRICAL APPARENT TIME TO EQUILIBRIUM, 
te, AND THE TIME TO EQUILIBRIUM, ts, 


PREDICTED BY THE KINEMATIC WAVE MODEL 


36 ts (sec) te (sec) 
0.05117 86.93 i-0.666 78.43 170-663 
0.0787 75yG7T 1-9-0060 72.15 470.687 
0.1045 68.46 i-0-666 66.29 i-0.695 
0.1583 59.22 170-666 60.29 i-0.702 


0.20791 54.12 i-0.666 52.08 i-0,684 


ct 


(8-8) 
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CHAPTER V 


STATE VARIABLE MODEL OF SURFACE RUNOFF 


5-1 Introduction 

Observations of runoff from natural watersheds and from laboratory 
catchments, as well as theoretical considerations of fluid motion support 
the view that rainfall and runoff are related by a nonlinear relationship. 
For this reason, accurate predictions of runoff resulting from wide 
variations of rainfall cannot be achieved by a linear model. It was shown 
in Chapter III, for the case of a single reservoir that a linear approxi- 
mation of a nonlinear process can be expected to be of an acceptable degree 
of accuracy only within a narrow range of fluctuations of the input vari- 
able (such as the rainfall intensity). There is a need to develop a non- 
linear model of the rainfall-runoff phenomenon based on relevant physical 
laws. 

The rainfall-runoff process is a part of a complex hydrological 
system which may have many inputs and many outputs, and these are related 
in a complicated manner. To analyze such a system, it is essential to 
reduce the complexity of the mathematical expressions, as well as to 
resort to computers for most of the tedious computations necessary in the 
analysis. The state-space approach to systems analysis, developed since 
around 1960, is best suited from this viewpoint (Ogata, 1970). The con- 
cept of state by itself is not new since it has been in existence for a 
long time in the field of classical dynamics (Ogata, 1967; MacFarlane, 


1970), and other fields. In hydrology, many investigators (Minshall, 1960; 
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Amorocho, 1967; Kirkly, 1967) recognized that catchment runoff is 
related not only to the rainfall input but also to the past rainfall- 
runoff history of the catchment. Both the antecedent index of precipi- 
tation and the base flow have been used as indexes of the state of the 
catchment for the purpose of predicting future runoffs. However, this 
practice has been largely based on intuitive considerations. 

In this chapter the state-space analysis is applied to the problem 
of simulating surface runoff due to rainfall distribution over an 
impervious uniformly sloping surface. The purpose of the analysis was 
to develop a mathematical model that could be used as a basis for a sub-— 


sequent empirical evaluation to simulate runoffs from natural watersheds. 


5-2 State Space Analysis of Systems 


The state-space method has been developed in modern system control 
theory. Modern control theory, states Ogata (1970), is contrasted with 
conventional theory in that the former is applicable to multiple-input- 
multiple-output systems, which may be linear or nonlinear, time-invariant 
or time-varying, while the latter is applicable only to linear time- 
invariant single-input-single-output systems. Also, modern control theory 
is essentially a time-domain approach and thus particularly suited for 
digital computer computations, while conventional control theory is a complex 
frequency-domain approach. 

From a mathematical viewpoint, the state-space approach to systems 
analysis and synthesis is the use of matrix and vector methods to handle 
the large number of variables which enter into such problems. As such, 
states Derusso (1967), these are not new methods, but rather they are the 


rediscovery of existing mathematical techniques. They aid considerably 
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in the solution of linear multivariable problems. More important, 
however, the state-space approach aids conceptual thinking about these 
problems, and nonlinear problems as well. Furthermore, it provides a 


unifying basis for consideration of linear and nonlinear problems. 


5-2-1 Dynamic Systems 


For the purpose of further discussion it is convenient to define 
dynamic and instantaneous systems. Schwartz (1965) states that a system 
is classified as instantaneous if its output at any instant t depends at 
most on the input value at the same instant but not on past or future 
values of the input. The input-output relation for an instantaneous 


system may be expressed by 


ze) t= Shu be lst) rk) 
where u = the input vector 

z = the output vector 

£ = the vector-valued function 


On the other hand, a system is said to be dynamic if the output 
at any instant depends not only on the present input, but also on at 
least some of the past input values. A system whose output at time t is 
completely determined by the input in the interval t-T to t, (T>0), is 
said to have a memory of length T. The input-output relation for a 
dynamic system may be expressed by 

Zeja= sc lel whist) (5-2) 
where x = the state vector 


the vector-valued function 


& 


All the relevant information concerning the past history of the system 


required to determine the response for any input is contained in the state 
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5-2-2 State and State Variables 

Ogata (1970) defines the state of a dynamic system as the smallest 
set of variables (called state variables) such that the knowledge of 
these variables at t = to, together with the input for t>t 9, completely 
determines the behaviour of the system for any time t>t ). Thus, the 
state of a dynamic system at time t is uniquely determined by the state 
at time to and the input for t>to, and it is independent of the state and 
input before to. 

If n state variables are needed to completely describe the 
behaviour of a given system, then these n state variables can be con- 
sidered to be the n components of a vector x(t) and accordingly such a 
vector is called a state vector. A state vector is thus a vector which 
determines uniquely the system state x(t) for any t>t,), once the input 


u(t) for t2to, is specified. 


Jvec3iState Space 


The n-dimensional space whose coordinate axes consist of the 
*], axis i cxgraxis, vicis xq axis'is cabled a*state space». Any state, 
determined by the state vector x(t), can be represented by a point in the 
state space. As the state of the system changes with time the locus of 
the points, representing the states of the system at different instants 


of time, follows a curve called the trajectory of the state vector. 


5-2-4 State Space Representation of Systems 


The state variable method of systems analysis applies to systems 


described by ordinary differential equations. Such systems are said to 
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be lumped. Systems which require the use of partial differential 
equations to describe their dynamic behaviour are said to be distributed. 
In the strict sense, states Schwartz (1965), all physical systems are 
distributed, but for practical purposes it is generally possible to 
approximate the behaviour of many by means of ordinary differential 
equations. 

The basic concept of the state variable method may be developed 


from the general equation 


(n) (n-1) ; 
Zetudy 2) een tt an=420h-an2 = u(t) (5-3) 
eda ve ee (n-1) e dn-l, (n)_ az 
where Z=—, Z ==, «-- 5 2 = ———- , = —— 
dt dt2 dtn-1 dt 


Equation (5-3) is general in the sense that it represents systems which 
may be classified in terms of the aj coefficients (i = 1, 2, ..., n) 
into the following three general categories (Langill, 1965): (1) linear 
Systems eit allsthescoeitictents, 1.6.4], a9, ..., ag, are constant, 
(2) linear time-varying systems, if any of the coefficients aj is a 
function of the independent variable time t, e.g. sint, 3t, e-2t, etc., 
(3) nonlinear systems, if one or more of the coefficients ay is a function 
of the dependent variable z, e.g. sin z, 3z, cos 2, etc. 

The state space representation of systems is based on the 
description of dynamical system behaviour in terms of n first-order 
differential equations. This may be done in a straightforward way, and 
without loss of generality, by adopting Moigno's device of introducing 
auxiliary variables (Moigno, 1844). To illustrate this technique on the 
example of Equation (5-3), let us define new variables x], x2, «+, Xp, 
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Xx] tA 
> dae 4 
(5-4) 
(n-1) 
Xn eZ, 


Since the specification of the values of these variables at a given 
reference time will determine a unique solution to the given equation, 
Equation (5-3), they are taken as the system state variables. Substitution 
of these variables then gives a set of n first-order differential 


equations which are equivalent to the given nth-order equation. These 


are 
Xx] = x9 
x2 = x3 
? (5-5) 
%n-1 = *n 
Xn SESank (Veen |X) ess ST kyr oO 
or 
Pee han (5-6) 
where 
X] 0 Ih 0 cae a O 0 
x9 0 0 1 v8 0 0 
P76] baal ce 1d eat id aig Re 
eet 0 0 0 ae ut 0 
Xn Pape Mane Agta Say i 


The output equation becomes 
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x] 
22> [110.4403 fs 
24 
or 
z= Cx (5-7) 
where 
C,=, [140 (es). 0) 


The first-order differential equation, Equation (5-6), is the state 
equation, and the algebraic equation, Equation (5-7), is the output 
equation. 

Whereas in the above example there was only one driving force, 
u(t), in a more general case m separate forces ul, u2, ..., Um will be 
present. In this more general case, the set of first-order differential 


equations will thus read 


X] 2 Ly C1, XD ss gee Uls Us wee 5 Um) 
RY = ERK Ky eee Mae UL oy es Um) 
a Feb) asics says UTS UD y os sig Un) 


Integrating Equation (5-8) from time t, to t gives 


£ 
xq(t)o= xp(te)at { feito soon plane xin tig Moelle eoadt, G=9) 
Co 


for i=1, 2, ..., n. From Equation (5-9) it can be concluded, states 
Elgerd (1967), that each of the n xj variables can be determined uniquely 
at any moment of time t if and only if: 


(1) Each variable x; is initially known, i.e. if n initial 
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conditions x;(to) are specified. 

(2) All m driving forces njnhes Lio? ort aten)mgare eee ee 
throughout the interval to to t. 

The n xj variables can thus be interpreted as carriers of the full 
information about the transient state of the system. Initially, at 
t = to, the total system state can be expressed by the n numbers x](to), 
X2(to), «+, Xn(to). Under the influence of the m driving forces and the 
bonds of interaction, given by the equations describing the system, the 
state of the system will change. The updated state expressed by the n 
numbers xj(t), x2(t), ..., xp(t) can be continuously obtained from 
Equation (5-9). In consequence of this interpretation, the state of the 
system is defined as the n-dimensional vector x(t), which has as its 


components the n numbers or state variables x1(t), x2(t), ..., xn(t). 


5-2-5 Second Order Systems and the Phase Plane 


The state space with only two dimensions reduces to a state plane, 
also known as a phase plane. The phase plane applies to those systems 
for which the differential equation is of the form 

“2 + az + aoz =u (5-10) 
Although the phase plane analysis is limited to first and second order 
systems it is important for several reasons. Firstly, it allows easy 
visualization, since a trajectory on a plane can be readily traced, and 
thus it enables better understanding of the concept of state and state 
space. Secondly, many practically important and dynamically complex 
systems can be reasonably well approximated by second-order systems. 
Finally, because of its simplicity, many significant results and 


conclusions have been obtained, and several important theorems concerning 
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system phase trajectories have been formulated (Hsu, 1968). 

The phase plane method is a method for obtaining graphically the 
solution of a first or second order equation given by Equation (5-10), 
which may be linear or nonlinear. The solution, i.e. the phase tra- 
jectory, is displayed in a phase plane which most commonly is the z-z 
plane. The procedure is basically the same as that for the state space. 
Equation (5-10) is rewritten in terms of a set of first-order differential 
equations by introducing the state variables x] = z, and x2 = Z. 


Equation (5-10) may be then written as 


x] = %2 
; (5-11) 
x9 =u-)}- ax = agXy 
or 
x9 2 dx» at u - a1)X9 - a9Xj (5-12) 


XL dxj x92 


Equation (5-12) defines the slope of the phase plane trajectory, dx2/dx,, 
in terms of functions x2 and x]. To find the curve which is the phase 
trajectory, this slope equation must be integrated. This integration may 
be performed analytically, graphically, or numerically. The trajectory 
does not show time information explicitly. If needed, however, the 
trajectory can be graduated in time units (Thaler, 1962; Ogata, 1970). 
For examples of the phase plane representation of second order 
systems the reader is refered to Appendix G. Definition of a conservative 
system is given and the graphical method of isoclines is explained. 
The analytical and graphical methods of solution are illustrated on 


the following solved problems: 
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Tl. Analytical Solution. 
(a) Determination of phase trajectories for a linear 
conservative system. 
2. Graphical Solution by the Method of Isoclines. 
(a) Determination of phase trajectories for a linear 
non-conservative system. 
(b) Determination of phase trajectories for a 
nonlinear system. 


Further examples can be found in Thaler (1962) and Ogata (1970). 


5-2-5-1 Phase Plane Analysis of Nonlinear Reservoir 


Approximation of outflow from a nonlinear reservoir by the trans- 
fer function of a linearized equation of the reservoir was presented in 
Chapter III. In the following the same reservoir is analyzed by the 
phene plane method. The dynamic behaviour of the reservoir is described 
by Equation (3-4), restated here for convenience 


ds a i 
7 KYS = i(t) (5-13) 


or more generally by Equation (3-3) 

Sa q(t) = i(t) 65210) 
where S is the storage, K is the reservoir constant, i is the inflow 
input, and q is the reservoir outflow. In the analysis of the given 
reservoir two separate cases may be distinguished depending on a zero 
or non-zero forcing function, i(t), during the considered time interval 
Ma oe i ie 


Case 1. Forcing function is zero, (iit] = 0). “The initial 
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89 
condition at time t = tp is given by S = So, and the equation of the 


reservoir is 


ds 
Sows =o (5-15) 


Taking S to be the state variable, the term dS/dt defines the slope of 
the trajectory in the S - t plane at a given instant of time t, for 
t>to. Putting N = dS/dt, where N is a constant, yields an equation of 
a family of isoclines 

S = (N/K)2 (5-16) 
The isoclines are seen to be straight horizontal lines. By choosing 
different values of N, such as Nj, N2, ..., Nx, a number of isoclines 
in the S-t plane can be determined. A particular trajectory can then be 
found by sketching a curve in the directional field of isoclines, 
starting from the point in the S-t plane given by the initial state Sp, 
at the time t = to. 

Three of the reservoir trajectories are shown in Figure 5-1. 

They fit very closely experimentally measured points. Since q = N (by 
comparison of Equations (5-13), (5-14), and (5-16)), both S and q can be 
determined from the same trajectory as functions of time. 

The method of isoclines offers an interesting possibility of 
"identifying" the system given by Equation (5-15) from a single record 
of outflow q. The discharge, q, from the reservoir is related to the 
storage, S, by Equation (3-2). 

q = K/S (5-17) 
Equation (5-17) defines a curve in the q-S plane. Slope N of this curve 


at any point P(q,S) is given by 
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(5-18) 


Suppose that the relation between q and S, given by Equation (5-17), is 
not known. Further suppose that a time record of outflow q is available 
in the range from q = qm to q = 0. Then for small increments of outflow 
q, Equation (5-18) may be approximated by 


y 2 Sg | Damen 


"aa ad 
AS or ge 
where qj, = the outflow at the time ty 
q2 = the outflow at the time ty 
AT = t2 - tj, is the time increment 


Equation (5-19) may be considered to define a family of isoclines 
in the q-S plane. Each isocline is a straight horizontal line given by 
the average outflow (q2 + q1)/2. When the isoclines are drawn the curve 
representing the relationship between q and S may be sketched and the 
parameters of the curve determined. 

This procedure was verified experimentally by measuring the dis- 
charge q from a laboratory setup described in Chapter III. Table 5-1 
shows the measured values of outflow q in discreet time intervals as well 
as the calculated average discharge qay and the slope N (given by 
Equation (5-19)) for.each time interval. The approximate curve deter- 
mined by the method outlined above, shown in Figure 5-2, fits very closely 
the exact curve plotted from Equation (5-17) in the same figure. The 
time, which appears as a parameter on the curve, can be evaluated as 
follows. The q-S relationship is approximated by straight-line segments. 
For a particular segment its slope N may be determined as well as the 


average discharge, qay, corresponding to the center of that segment. 
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TABLE 5-1 


PHASE PLANE ANALYSIS OF A LABORATORY RESERVOIR FOR THE CASE OF 


Time 
(sec) 


0.00 

5.95 

7.86 
11.7 
1 Bs i W 
19.65 
25008 
PA Boe RI 
31.44 
Sheep 
$9230 
43.23 
47.16 
51.09 
55.02 
58495 
62.88 
66.81 
70.74 
74.67 
78.60 
82.53 
86.46 
90.39 
94.32 
98525 
102218 
106.11 
110.04 
£13297 
117,290 
E2183 
125.76 
129'.69 


RECEDING OUTFLOW AND ZERO INFLOW: 


ISOCLINES CALCULATED BY EQUATION (5-19) FROM 


MEASURED VALUES OF DISCHARGE 


Discharge 
(cm3/sec) 


47.7 
46.3 
45.0 
43.7 
42.4 
41.1 
Be 
38.7 
37.4 
BOe2 
34.9 
6 be 
32) 
3Ee3 
30.50 
28.8 
27.2 
Ole 
24.9 
23.7 


Average Discharge 
(cm3/sec) 


47.0 
45.7 
44.4 
43.0 
41.8 
40.5 
B9e0 
38 oL 
36.8 
350 
34.3 
Shap 
SNe, 
30.7 
29.4 
ZOeL 
269 
20 
24.3 
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Figure 5-2. Phase plane solution of storage equation 


of a nonlinear reservoir for the case of zero inflow. 
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The displacement Aq = qo - qj can be measured over the length of the 
segment. Then from Equation (5-19) the time increment, AT, for that 


segment is 


A 
Ndav 


Time for any point j on the phase trajectory may be calculated by 


(5-20) 


summation of the time increments, AT;, starting from the intial point, 
J 
2 ES | 2 aa Saag 5 


Equation (5-19) may be rearranged into 


2-NAT 
q2 = 41 DeNAT sea 


which has the form of a well known depletion equation defining a 
recession curve of a hydrograph, usually written as 

dt = dt-1 K& (5-22) 
where Ky is a recession constant. Equation (5-22) is often used in 
hydrograph analysis. It can be seen, by comparing Equations (5-21) and 
(5-22) that Ky can be constant only if N is a constant. However, a constant 
N implies a linear relationship between q and S (by Equation (5-18)) 
which is a questionable assumption with respect to watershed outflow. 
Thus, Equation (5-22) must be regarded as a gross linear approximation. 

Case 2. Forcing function is not zero (i(t) # 0). This case was 
solved in Chapter III by the transfer function method for step inputs, 
using a linearized model of the reservoir. The initial conditions at 
time t = to were: S = Sy, ix = qx, and for any time t>t > the input was 
i=i, +I. The asterisk signifies a steady state value and I is the 


step increment of the input. Using the phase plane analysis the same 
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problem may be solved graphically as follows. 
The system equation is given by Equation (5-13), i.e. 


= = i(t) - K’S (5-23) 


The phase trajectory in the S-t plane may be approximated by short 
straight-line segments tangent to the trajectory. The slope of each 
tangent line segment is N = &, which can be calculated successively 

for small time (or storage) increments from Equation (5-23). This process, 
in principal similar to numerical solution of a differential equation, 

is organized in the following manner. Starting from the initial con- 
ditions at time t = to the trajectory slope N at this time is calculated 
by Equation (5-23). A tangent line segment, having the slope N, is 
drawn. The length of the segment must be estimated so as to approximate 
closely the trajectory. A new value of S is read at the end of the 
segment. Substitution of this new value of S into Equation (5-23) yields 
a new value of the slope, No. Another tangent line segment is drawn, 

and the whole process is repeated. The tangent line segments form an 
envelope to the trajectory. 

Responses of the reservoir to two step function inputs (the same 
problem as described in Chapter III) determined by the above method are 
compared to experimentally measured responses in Figure 5-3. The 
agreement between calculated and measured trajectories is very good. 
Further comparison of results obtained by the phase plane method and 
those obtained by the transfer function of a linearized reservoir 


indicates clearly the superiority of the former method. 
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5-3 State Variable Model of Spatially Varied Overland Flow 


Overland flow on a watershed is essentially a three-dimensional 
free-surface unsteady spatially varied flow phenomenon. However, unless 
there is a well defined channel, surface runoff can be approximated as 
two-dimensional (Yen, 1968). Thus, for the analysis of overland flow 
resulting from rainfall, consider an idealized impermeable surface shown 
in Figure 5-4. The surface is of length L, slope Sg = sin 6, and of 
unit width. The x axis is taken along AB with the origin at A. When 
rain of a constant intensity i begins to fall the discharge changes with 
time and the flow is unsteady. A time of equilibrium will eventually be 
reached. At this time the discharge is equal to the rate of rainfall 
and the runoff becomes steady. When rainfall ceases, the subsiding 
runoff becomes unsteady again. The objective was to develop a 
mathematical model (using the state space formulation) of this one- 


dimensional, spatially varied unsteady surface flow. 


5-3-1 Unsteady Spatially Varied Flow - Basic Equations 


The equations of two-dimensional unsteady spatially varied flow 


are the continuity equation 


oy, Saye i cos 6 (5-24) 
ot ox 


and the dynamic equation 


oU 4 ota ey 2 ie we cif co 
rte relia ai g(So-S¢) ae (5-25) 


In the above equations x is the longitudinal direction of the flow along 
the x-axis of the runoff plane. The average longitudinal velocity 
component of the flow is U. The lateral inflow is i, in terms of volume 


per unit time t per unit surface area on the runoff plane. This lateral 
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Figure 5-4. Definition sketch for surface flow 


over a sloping plane. 
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inflow has a horizontal velocity component u in the x-direction, y is 
the depth of flow normal to the plane and the slope of the runoff plane 
is So, S¢ is the friction slope. 

These equations were first derived by de Saint-Venant in 1871, 
and therefore, are sometimes referred to as the de Saint-Venant equations. 
Derivation of the equations with their assumptions may be found elsewhere 


(Rouse, 1957; Chow, 1959; Chen, 1968; Yen, 1968). 


5-3-2 Unsteady Spatially Varied Flow - Approximation by Lumped Parameters 


While the de Saint-Venant equations describe the unsteady spatially 
varied flow in a rigorous manner, their practical usefulness for 
simulating overland flow from watersheds is somewhat limited. This is 
because of the complexity of the equations. Their analytical solutions 
have been obtained only for restricted, simplified cases. The only 
dvertap i general method of solution are finite difference techniques 
(Chow, 1959; Strelkoff, 1970; Brutsaert, 1971). However, these methods 
have a major disadvantage for continuous simulation since substantial 
amounts of computer time are needed. Furthermore, areal variations in 
the rainfall inputs as well as the surface parameters of natural water- 
sheds, such as the infiltration rate, surface roughness, length and 
slope, require the use of average values. Thus, while reasonably 
accurate methods are justified, the accuracy to be gained by using 
finite difference methods for overland flow, states Crawford (1966), is 
still subject to question because of the limited accuracy of the basic 
data. 

It is thus desirable to develop a model based on simplified 


mathematical relations. This can be achieved by approximating the 
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distributed system by lumped parameters. In the following the con- 
tinuity and dynamic equations of surface runoff, based on lumped para- 


meters, are derived. 


5-5-2-1 the Continuity Equation 


Equation of continuity of overland flow, considered as a 


distributed system shown in Figure 5-5(a), may be written as 


Od 4. 2Y as = 
yee am 1 cos (0 (5-26) 


where q = Uy is the discharge per unit width. Dynamics of this distri- 
buted system can be approximated by a lumped parameter model, shown in 
Figure 5-5(b), in which the runoff plane is substituted by a series of 
interacting reservoirs. The distributed parameters i(x,t) and q(x,t) 
indicating the rate of rainfall and discharge respectively, have been 
lumped into a finite number of input and output variables, Pj and 04 
respectively, where j = 1, 2, ..., n. This lumped model is of a 
"distributized" character since it consists of a number of lumped para- 
meter units or subsystems. The equation of continuity, generally valid 
for each of the reservoirs, can be written in terms of the ordinary 


differential equation 


dSj _ dhy _ a 

qe 7 Xicos 8 As =P, = (0j = 04-1) (5-27) 
where Sj = the storage in the j-th reservoir per unit width 

Xj = the length of the j-th reservoir in the x-direction 

hy = the depth of water in the j-th reservoir 


The lumped input Pj is given by 


Ps =i Xj cos 6 (5-28) 
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and the output 0; is to be determined by the dynamic equation. Continuity 
of the whole model consisting of n reservoirs is then described by the 
set of n ordinary differential equations 


EK cones Seite 
era: WEA Lee) 


dh 
Xo cost P9 = (O59 = 04) 


: (5-29) 


dh 
Xn cost=— = Ph - (On - On-1) 


More refined approximation to the distributed system may be 
achieved by routing the flow through n hypothetical reaches instead of 
reservoirs, as shown in Figure 5-5(c). Mechanism of outflow from these 
reaches is assumed to be such that the water surface remains parallel 
to the bottom at all times. The set of continuity equations 


characterizing the model is 


dy1 

X] ate = Py - Oj 
d 

X9 22 P20) 104) 

e ° e (5-30) 
dyn 

Xp a Pn = iva On-¥2 


The approximation to the distributed system becomes progressively 
better as the number of reaches, n, in the lumped model increases. This 
can be proved by letting n increase to infinity, i.e. Xj approaches zero 


(j = 1, 2, ..., n). The left hand side of any of the equations (5-30) 


lim (X; 21) = 6x°¥. (5-31) 
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Figure 5-5. Surface runoff over a sloping plane: (a) definition 
sketch, (b) simulation by a 4th-order lumped system of interacting 
reservoirs in series, (c) simulation by a 4th-order lumped system 


consisting of hypothetical uniform flow reaches. 
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and the right hand side is 
: : 3 
lim [P. - (0; - 05-1)] = idx cosé - 5i6x (5-32) 
X70 
J 
By comparison of both sides one obtains 


cb Dea: Cy a 
Wes we (5-33) 


which is the continuity equation, Equation (5-26) of the distributed 


system. 


5-3-2-2 The Dynamic Equation 


Since the lumped model assumes uniform flow within every reach, 
discharge may be calculated by any type of resistance formula. In 
general terms the dynamic equation has the form 

O = aym (5-34) 
The parameters a and m are given as follows: 


(1) Laminar flow (derivation can be found in Chow (1959)): 


S 
a = BS, 


(2) Turbulent flow: 


m= 3 


Chezy equation: a = C182, m = 3/2 


Manning equation: a = CoS’, m = 5/3 


where g = the acceleration of gravity 
S = the total energy line gradient 
v = the kinematic viscosity of water 


Cj = Chezy resistance coefficient 


C2 = Manning resistance coefficient (1.49/n) 


The Reynolds criterion for stability of flow on changing from laminar to 


turbulent may be used to determine the appropriate formula. More 
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discussion about the Reynolds criterion follows in Section 4. 


5-3-2-3 The Energy Equation as a Bond of Interaction 


Definition of the slope S in the parameter a is of basic impor- 


tance. Consider a prismatic channel in which the flow is gradually 


varied, shown in Figure 5-6. The total energy at the channel section 
1 ds 
” 
Hy = YI 
i =.2) + y] cosd oF (5-35) 
& 


and at the channel section 2 is 


U 
H2 = z9 + y2 cosd + a2 (5-36) 
2g 
The energy gradient, denoted by Sf, between the sections 1 and 2 is 
ie) 
H7-H Fie: 
Sf = ee = an ®Y cose 7s 6. (5-37) 
The water surface slope, Sw, is 
Sw = Az 4 Ay cosé (5-38) 


L L 
and the channel bottom slope, So, is defined as 


So = = sind (5-39) 


In uniform flow, S¢ = Sy = So = sind (Chow, 1959). 

Now, if the slope S, is used in the dynamic equation, Equation 
(5-34) in conjunction with the set of continuity equations, Equation 
(5-30), to simulate surface runoff, the following situation will result. 
With reference to Figure 5-7(a) the actual steady state flow at cross- 
section j is 

‘ m 
qd; = Li= ay; (5-40) 


where L is the distance in the x-axis direction from the origin to the 
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Figure 5-6, Energy in gradually varied overland flow. 
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Figure 5-/. Definition sketch for comparison of: (a) actual 


flow, and (b) flow simulated by noninteracting reaches. 
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j-th cross-section. The same discharge will appear in the simulated 
flow shown in Figure 5-7 (b), assuming that the dynamic equation is 
correct. Thus, the depen of the actual and simulated flow at the j-th 
cross-section will be the same, i.e., qj = Oj, and yj = ¥4 The 
geometrical interpretation of this result is that the simulated steady 
State profile will approximate the actual profile by a step-like curve. 
The depth at the lower boundary of every reach will correspond to the 
actual depth but because greater storage is required in the model to 
reach the steady state, the time to equilibrium will be longer for the 
simulated flow than for the actual runoff. 

While the use of the channel bottom slope, So, in the model 
simplifies greatly the computation, it has one serious disadvantage. 
The model is not capable of simulating backwater curves because there is 
no interaction between the reaches. This can be corrected by using the 
energy gradient, S¢, in the dynamic equation. Definition of this 
slope for the lumped parameter model is apparent from Figure 5-6 


Considering the cross-section A, the slope of the energy gradient is 


Zoree 
Ay eee cos ad U}-U 
1-22 et (y1l-y2) fi 1582 


L L 2eL atl) 


Equation (5-41) may be considered as a bond of interaction between the 
two reaches. Similar equations apply to all interfaces between the 
reaches of the model. It is conceivable that in some cases of gradually 
varied flow the last term in Equation (5-41), representing the difference 
in velocity heads, may be negligible and the slope of water surface 


Z1-22 , cos§ (y1-y2) 


: (5-42) 


Sw = 


can be used in the dynamic equation. 
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5-3-3 The State Variable Model 
The lumped parameter model of surface runoff based on the con- 
tinuity, dynamic, and energy equations given by Equations (5-30), (5-34), 
and (5-37), may be described mathematically in a compact matrix form 
using the state space representation of systems. The state and output 
equations of the lumped parameter surface runoff model are defined as 


follows. The state equation is 


y=Ay+BP Co 
where 
yi 
al mae 0) 0 eee) 0 
yhed oy 
x2 eee sae 
Oa) Se 0 0 0 
y = 3 A = e 
a j 
n— 
Yn-1 0 0 Ries Rey § vee bits eres 0 
ye-1 ya 
0 0 0 BS Sas 
Yn Xp-1 n 
I 
x] 0 0 ‘ 0 0 
1 
0 Te 0 0 0 
B = 
uh 
0 0 0 5 ome |!) 
Xp-1 
1 
0 0 0 0 X, 


and P is the precipitation input vector. 
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The output equation is 


Ore Sow (5-44) 
where 
2 

Oj ayy 0 0 ° e ° 0 0 
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On-1 0 0 Oe etl aa) eae Re voL aL. MO 
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Operation of the lumped parameter runoff model, given by Equation 
(5-43) and (5-44) may be visualized in the following steps: 

(1) The runoff plane is subdivided into a number of n sub-areas 
which constitute the elements of the lumped model. The flow within each 
element is assumed to be uniform at all times. 

(2) The outflow from each of the n elements of the model is 
determined for small increments of time, AT, by solving simultaneously 
the set of continuity equations, Equation (5-43), followed by 
simultaneous solution of the set of output equations, Equation (5-44). 

(3) The input into each element j (j =1, 2, ..., n) at the 
beginning of a particular time interval AT is assumed to be an impulse 
of a magnitude equal to the sum of the rainfall input, Pj, the surface 
inflow, 0j-], and the surface outflow, 0;. The impulse input causes 
an instantaneous change in the depth of flow. However, the water 
surface remains parallel to the bottom. Thus, the flow within each 


element may be characterized as uniform unsteady flow. 
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(4) The response of each element to the impulse input determines 
the outputs and the state of each element at the end of the time 
interval AT; the new states are the initial conditions for the following 
time increment. 

(5) Each of the elements operates essentially as the nonlinear 
reservoir analyzed in Chapter V, Section 2-5-1, by the phase plane 
method. Since the lumped parameter runoff model consists of n elements 
the solution is performed in n-dimensional state space rather than in a 


two-dimensional phase plane. 


5-4 Modification of the Dynamic Equation for Raindrop Impact Effect 


In overland flow with rainfall, drops after striking the surface 
water layer generate a local disturbance, which results in high resis- 
tance to flow. Although the raindrop impact effect has been observed, 
its mechanics as yet has not been fully explained. The problem of 
rainfall impact effect, as related to surface flow, was investigated 
in the present study for the purpose of testing the state variable 


runoff model. 


5-4-1 General 

The retarding influence of rainfall impact effect on surface flow 
has been observed by many investigators. Early experimental studies by 
Izzard (1943; 1944, 1947) confirmed the existence of the effect of drop impact, 
which demonstrated itself by a momentary increase in the discharge after 
cessation of rainfall. This phenomenon was attributed to the fact that 
during rainfall the amount of detention is greater, due to the raindrop 


impact effect, than the amount required to cause the same rate of runoff 
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BEL 
after cessation of rainfall. Wenzel and Yoon (1971) measured surface 
flow depth with rainfall and recognized that it is substantially deeper 
than that without rainfall at the same flow rate. 
The rainfall impact effect also influences the transition of flow 
from laminar to turbulent region. Reynolds number, which is used as a 
criterion for distinguishing the laminar and turbulent flow, is defined 


for overland flow as 


Re = aa (5-45) 
where U is the average flow velocity. Yu and McNown (1964) found from their 
analysis of the Los Angeles airfield drainage data that the transitional 


range of the Reynolds number is 200 and 1000 for flow on a concrete 


surface with and without rainfall respectively. 


5-4-2 The Relationship Between Rainfall Intensity and Friction Factor 


The friction factor is a measure of the resistance to the flow, 
and for flow in pipes or open channels it gives the relationship between 
the energy gradient slope (S¢), mean velocity (U), and the depth (y) or 
hydraulic radius (R) of the flow. It can be expressed in a number of 


different forms including the Darcy-Weisbach equation 


8eyS 
pcs SY°f 
u2 


(5-46) 
For flows in pipes and in rigid boundary open channels the friction 
factor depends primarily on the relative roughness (e/y) in the 
hydraulically rough regime, on the Reynolds number of the flow (Uy/v) in 
the hydraulically smooth regime, and on both these parameters when the 


flow is transitional between the two regimes. However, it appears that 


in case of sheet flow under rainfall the friction factor also depends to 
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Pe 
some extent on the rainfall characteristics such as the rainfall 
intensity, size and spacing of raindrops, and their impact velocity. 

Woo and Brater (1962) carried out sheet flow experiments with 
controlled rainfall in which they studied the relationship of the friction 
factor of flow under rainfall, f;, Reynolds number, Re, and the runoff 
surface slope, So. They made the following conclusions: 

(1) The raindrop impact effect is greater for small depths and 
small values of Re. 

(2) For rough surfaces the relationship between fy and Re differs 
from the steady flow relationship more than for smooth surfaces (because 
turbulence is more easily generated on rough surfaces). 

(3) The fy values for a particular Re decrease with increasing 


< 0.006), but do not follow this rule for 


slope for small slopes (So 
slopes steeper than 0.008. 
(4) For small slopes, values of f,; are larger than the corresponding 
values of f for uniform flow without rainfall. 
Yoon and Wenzel (1971) confirmed the findings of Woo and Brater 
in their report on mechanics of sheet flow under simulated rainfall. 
They found that for Re <2000 the rainfall intensity increased the 
resistance of flow at a constant Reynolds number. Above R, = 2000, the 
intensity effect decreased rapidly until it became insignificant. For 


R: 


a < 1000 the results could be described for practical purposes by 


straight lines, parallel to the theoretical f = 24/Re which is valid for 
laminar flow. The parallel lines were thus given by a general relation- 
ship £ = C/Re, in which the authors expressed C as a function of channel 


slope and rainfall intensity. 
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The same authors also studied effects of changes in the impact 
velocity, from 8.8 ft/sec to 22.3 ft/sec, and raindrop spacing of one 
and two inches, on the friction factor. These changes were found to be 
insignificant. Similarly, the Froude number did not yield any significant 


correlation with the friction factor. 


5-4-2-1 Dimensional Analysis 


The objective of dimensional analysis is to express the behaviour 
of a phenomenon in the form of functional equations relating non-dimen- 
sional parameters. These equations can then be evaluated empirically 
to determine the relative importance of individual parameters and to 
formulate the relationship between these parameters. 

The method of dimensional analysis may be used in the development 
of the relation defining the friction factor, fr, of flow under rainfall. 
Considering the case of runoff due to rainfall from an impervious surface, 
the physical variables which most likely effect the flow phenomenon may 
be expressed in a functional form as 

Be = Fa (Us Vp Sos s0s8>elsd.Uren A) (5-47) 
where f, = the friction factor of flow under rainfall 

U = the local mean flow velocity 

y = the local depth measured normal to the channel bottom 


So = the channel slope or sind 


6 = the angle between the x direction and the horizontal 
u = the dynamic viscosity of water 

0 = the density of water 

g = the acceleration of gravity 


¢€ = the equivalent surface roughness 
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i = the rainfall intensity 
d = the mean drop size 
Uy = the raindrop impact velocity 
n the dimensionless parameters describing the drop impact 
r eleeen and shape 

The last four terms in Equation (5-47) were neglected in further 
analysis because they were kept constant during experiments, and further- 
more were shown by Yoon and Wenzel (1971) to have little significance. 
Thus, the functional relationship of Equation (5-47) simplifies into 

foe Ena (U.S, el Sine 5 1) (5-48) 
Application of Buckingham theorem of dimensional analysis (Blench, 1969) 
leads to 


frerirn) (CUZ. ia eel es (5-49) 
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The effects of the Froude number and the slope So were considered to 


have only secondary effects, thus giving 


fp=m@,£, 5 (5-50) 


y * Udy 
Equation (5-50) is a general functional relationship describing the flow 
from an impervious surface caused by rainfall. Since all flows in the 
experimental part of the present study occurred within the laminar 


region, further attention was given specifically to the relationship 


5-51 
fr = Fn &, am ( ) 
y 


which applies to laminar flow under rainfall 
For laminar flow without rainfall the functional relationship is 


known to be (Blench, 1969) 
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f= Fn (2) (5-52) 


which can be expressed by a general equation 


f== (5253) 


A theoretical value of C can be found from the Darcy-Weisbach equation, 
written for an infinitely wide channel 


8gSoy 
2 


(5-54) 


which is valid for both uniform and nearly uniform flows (Chow, 1959). 
Comparison of Equation (5-53) and Equation (5-54) gives 


8gSoy" 
ieee csr eee 


The theoretical mean velocity of laminar sheet flow (Chow, 1959) is 


gSoy 
ge (5-56) 
which after substitution into Equation (5-55) gives 
C = 24 (5-57) 


Solving Equation (5-55) for the mean velocity and substituting into the 
relation for specific discharge, 0 = Uy, yields 


not 8gSoy3 


Cy (5-58) 


It was hypothesized that the raindrop impact effect on laminar 
flow can be included in the friction factor f;, expressed as 


eS 


fr ce Re (5-59) 
where 
S.v2 4 
Gis Ceeay =” Fn (2 (5-60) 
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or 
Cr = C+ Cy (5-61) 


Thus, Equation (5-58) is modified for the case of runoff under rainfall 


into 
8eSoy3 
O = BELO (5-62) 
Crv 
The rainfall intensity term 
iv 


becomes zero when rain ceases, and Cy is then equal to C. The impact 
effect diminishes with increasing depth and velocity of flow, whereas 
increasing rain intensity and higher viscosity of the fluid tend to 
enhance the impact effect. These relations are in accordance with 


intuitive expectations. 


5-4-2-2 Empirical Evaluation 

The functional relationship given by Equation (5-63) has to be 
evaluated empirically. Because of very small depth of flow involved in 
the experiments, the depth and velocity could not be measured directly. 
An approximate method was therefore used, which was based on measured 
amounts of water detained on the runoff plane during steady flows, for 
various rainfall intensities. 

With reference to Figure 5-8 a steady state discharge Q from 


a runoff plane of length X and width B is 


Ora Bilas (5-64) 
and also 
3 2So 
Q = Ba Pee » where gq = oe (5-65) 
Cr Crv 
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Figure 5-8. Definition.sketch for steady discharge condition 


(Q = BiX = Bay?/Cr) from a plane of length X and width B. 
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Solving both Equation (5-64) and Equation (5-65) for y gives 


Te (5-66) 


Y= 


which is an equation of the water surface profile at steady state flow 
conditions. 
The volume of water temporarily detained on the runoff plane is 


given by 


X 
S-= B | ydx (5-67) 
0 


Substituting for y from Equation (5-66) yields 
X 
eS Cou: xt/3 ax (5-68) 
0 
Equation (5-68) is an integral equation expressing the amount of detained 
water, S, during the steady runoff as a function of distance x, rainfall 
intensity i, and the coefficient Cr. The coefficient Cy depends on the 
rainfall intensity i and the depth of flow y, and thus it is an implicit 
function of distance x. Exact integration of Equation (5-68) is therefore 
not possible. 
To evaluate the coefficient C, approximately, an "average" depth 
of flow, y, in the reach of length X was substituted into Equation (5-67), 


Solution of Equation (5-65) yields an equivalent coefficient Cy 


id 25 
Cr = Bo (5-69) 


The equivalent coefficient Cy is representative of an equivalent uniform 


steady flow of depth y and discharge Q within the reach of length X. The 
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method thus involves approximation of actual non-uniform flow by an 
"equivalent" uniform flow, with the same output discharge, and one might 
classify the method as lumping of parameters. 

Experimental evaluation of the parameters ee Sik and y for 
different rainfall intensities and for three slopes of the runoff plane 
was carried out by dividing the runoff plane into four reaches of equal 
length X = 19.6 cm, as shown in Figure 5-9. In the Fare measurement 
of discharge Q] (Figure 5-9(a)) the rain was allowed to fall on the 
first reach only. The amount of detained water, S], at steady flow was 
determined graphically from the vertical distance between the recorded 
curves of cumulative input and output, after the steady flow was reached 
as illustrated in Figure 5-10. 

In the second measurement the rain was allowed to fall on the 
first and second reach, and the total amount of stored water, Sy, was 
determined by the same graphical method. The term S9, representing the 
amount of stored water on the second reach, was found by subtracting 
Sj from St. By successive adding one more reach in each measurement, 
the amount of stored water in each reach was determined from the general 
formula 

S, © Stu= Saad (5-70) 

Measurements were repeated for several rainfall intensities, and 
for the three slopes: 12° (S, = 0.20791), 6° (S, = 0.1051), 2°56' 

(Sg = 0.05143). The intensity factor, Cy, was evaluated from the 
rearranged Equation (5-61) 
Cy = Cy — 24 (5-71) 


The Cy factor was plotted in Figure 5-11 against the dimensionless 
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Figure 5-9. Definition sketch for determination of storage 
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Figure 5-10. Determination of the amount of stored water, ST, 


on the runoff plane by means of mass curves of inflow and outflow. 
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Figure 5-11. Best fit line of the relationship between the 


intensity factor Cj and the dimensionless parameter iv/Uty. 
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Figure 5-12. Friction factor, fy, for laminar sheet flow 


under rainfall versus Reynold's number Re. 
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parameter » Equation (5-63), on a semilogarithmic paper with Cj 


uly: 
Rev 
values on the arithmetic scale. A best fit straight line, fitted to the 


data by the least square method, is given by the equation 


oe. LO Fad) 
Go = ln 5-72 
Epsane 17h ae ae 


with the correlation coefficient r = 0.922. Substituting the above 
expression into the equation for friction coefficient under rainfall, 


Equation (5-59), yields 


Uris Be8oy- | 10> 
$ExT Rel Uv ¥ Teese! Rev! (5-73) 


or for a laminar sheet flow 


1 1 10? i 
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A summary of all experimentally determined, and calculated para- 
meters is given in Appendix C. Values of fy and Re, plotted in the form 


of a friction diagram, are shown in Figure 5-12. The plotted points, 


24 


following approximately a parallel line to the theoretical f = . 
e 


’ 
indicate higher values of the friction coefficient then for laminar flow 
without rainfall. The analyses presented in this thesis is based on the 
premise that laminar flow prevailed on the model independent of rainfall 
impact effects and therefore the equation: Qay? was valid. All values 

of Re based on depths of flow calculated from Equation (5-67) (or y =S/BX) 
and the average velocities of flow calculated from U = Q/By were < 30. 
This value of Re is well within the laminar flow limits suggested by other 
investigators (Yoon, 1971; Schreiber, 1972). Further, the results pre- 
sented in Figure 5-12 give additional credence to this assumption because 


of the close agreement between the experimental results and the 


theoretical curve: f = 24/Re. 
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5-4-3 Comparison of Results with Existing Knowledge 


Examining Equation (5-74), one might expect that for the same 
Reynold's number the friction coefficient should vary with the rainfall 
intensity. This was found to be the case by Yoon and Wenzel(1971) for 
Re < 2000. The authors presented a friction diagram, where experimentally 
found values for fy followed parallel lines to f = a, each parallel 
for a given rainfall intensity. In the present study such distinction 
due to different intensities could not be made, most likely because of 
the approximate method used to determine fy and Re, and because the 
small depth of flow, ranging from about 0.2 mm to 0.5 mm. It is 
probable that for such small depths the effect of varying intensity of 
rainfall is not apparent. 

The experiments conducted by Yoon and Wenzel also showed, that 
for Re>2000, the rainfall intensity effect decreases rapidly with 
increasing Re until it becomes insignificant. The authors suggested 
that the retardation effect of the rainfall intensity on flow can be 
explained on the basis of momentum considerations. If raindrops strike 
almost normal to the water surface, the component of momentum carried 
by raindrops in the mean flow direction is almost negligible. Therefore, 
a portion of mean flow momentum has to be transferred to accelerate the 
raindrop mass from almost zero velocity in the x-direction to the point 
velocities at various depths. The velocity of mean flow is thus some- 
what retarded below the velocity that it would have if no rainfall were 
present. Because the raindrops strike the water surface first and then 


diffuse vertically as well as laterally, transfer of mean flow momentum 
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is greatest near the surface at high Reynold's numbers than for the low 
ones and the fact that the rate of raindrop impact is constant for a 
given intensity and drop size, the relative magnitude of transfer of 
mean flow momentum decreases with increasing Reynold's number. 

The findings in the present study are in general agreement with 
results of recent investigation conducted by Schreiber and Bender (1972), 
who obtained overland hydrographs from a laboratory plane catchment 
8 ft by 16 ft, with a lightweight-concrete surface. The rainfall was 
generated by rows of nozzles. The hydrographs were simulated by a 
mathematical model based upon the simultaneous solution of the continuity 
equation, Equation (5-26), and the dynamic equation, Equation (5-34), by 
a finite-difference method. The resistance coefficient K (corresponding 
to the coefficient Cr in the present study) in Equation (5-59) was 
determined for each hydrograph by parameter optimization, i.e., the 
unknown parameter K was varied systematically to obtain the best fit 
between the computed and observed hydrographs. However, this procedure 
assumed that K is constant over the entire flow plane at a particular 
instant in time. 

The coefficient K was plotted against various dimensionless 
parameters. The most successful plottings involved K, Weber number W, 
and precipitation number P. The latter parameter is the volume rate of 
precipitation divided by the product of mean velocity of flow and the 


square of the depth of flow. Values of K appeared to be dominated by P 
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during rainfall and by W after rainfall. 

The magnitude of K ranged from 29 to 40 for the rainfall 
intensity around 3 in/hr, and from 20 to 25 for the intensity around 
l in/hr. Optimal values of K for the recession periods ranged from 13 
to 18. The resistance coefficient, C;, proposed in the present study 
ranged from 24 to 32 for the rainfall intensity range from 0 to around 
4.5 in/hr. The value of Cy for the recession periods was 24. The 
coefficient Cy is based on physical reasoning whereas the parameter K 
was determined more or less by a trial method. The most significant 
objection to the method discussed above is the assumption of constant K 
over the entire flow plane at an instant in time. In fact, as the 
authors pointed out themselves, the friction factor f, and therefore 
also K, varies with Reynolds number, and thus K cannot be constant along 
the flow plane for non-uniform flow. 

The dimensionless precipitation number, P, which the authors 
correlated with the coefficient K does not include the kinematic 
viscosity of fluid. The kinematic viscosity was considered significant 
in the present study. The parameter P also includes a rather awkward 
quantity - the volume rate of precipitation. The expression for the 
friction coefficient derived in the present study, Equation (5-73), is 


considered more meaningful physically. 


5-5 Comparison of Predicted and Measured Hydrographs 


The detailed description of experimental techniques was given 
in Chapter II. Basically, the laboratory catchment hydrographs were 
recorded in the form of discrete time series, which were further pro- 


cessed by a computer. The end product of the computer analysis were plots 
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of runoff hydrographs delineated by a series of points plotted at the 
sampling intervals of the discrete data. Each runoff event for a 
particular set of conditions (rainfall intensity, temporal and spatial 
rainfall distribution, slope of the runoff surface, and initial wetness 
of the surface) was repeated seven times to obtain a better estimate 
of the measured values. 

Recorded hydrographs were simulated by a seven order state variable 
runoff model (thus the length of each of the seven reaches was 13.06 cm 
(5.14 in)). The dynamic, or discharge formula of the model was given 
by Equation (5-34), which included the derived relationship for raindrop 
impact effect on laminar flow, Equation (5-74), in the form of an extended 
friction factor, fy = Cr/Re. Figure 5-13 shows the effect the factor 
fr has on the solution as compared to the solution when the raindrop 
impact effect is not included, i.e., for f = 24/Re. The theoretical 
value of the resistance coefficient for laminar flow with no rainfall, 

C = 24, was used for simulation of all recession curves. 

The effect of a number of stages (reaches) in the mathematical 
model (i.e., the order of the simulating system) on the model's response 
is illustrated in Figure 5-14. As the number of stages increases, the 
S-shaped response of the model (series of interacting nonlinear stages 
with a distributed input) steepens and approaches a limit, given by the 
response curve of a completely distributed system, as was proven in 
Chapter V, Section 3-2-1. 

The experiments were divided, according to temporal and spatial 
distribution of the rainfall input, into the following three categories: 


(1) Experiments with uniformly distributed rainfall input. 
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(2) Experiments with temporally varied rainfall input. 

(3) Experiments with spatially varied rainfall input. 
Plotting of the experimental data and the state variable model simulation 
into a common plot for each event has been used as a mee of displaying 
and appraising the results. Statistical analysis of the experimental 


data for accuracy are given in Chapter VI. 


5-5-1 Experiments with Uniformly Distributed Rainfall Input 


Hydrographs of the surface runoff, generated by temporally and 
spatially uniform rainfall inputs, formed the bulk of the laboratory 
catchment experiments. Three factors were varied in the experiments: 
(1) slope of the runoff plane (approximately 3, 6 and 12 degrees), 

(2) rainfall intensity (low, i = 40.0 mm/hr; medium, i = 75.0 mm/hr; 

high, i = 115.0 mm/hr), and (3) duration of the rainfall input (two 

short durations, 8 to 31 sec, which were shorter than time to equilibriun, 
and one long duration of 200 sec). 

Exact rainfall intensity in every event was calculated from 
steady flow conditions resulting from the long duration (200 sec) rain- 
falls. 

The laboratory catchment experiments with uniformly distributed 
rainfall inputs were carried out in the following manner. For each of 
the three slopes the catchment hydrographs were measured for three 
rainfall intensities (low, medium and high). Further, hydrographs 
resulting from three rainfall durations were measured for each of the 
rainfall intensities. Thus, 27 different runoff events existed, and 
since each was repeated seven times, the total number of recorded hydro- 


graphs was 189. Summary of the runoff events is given in Table 5-2. 
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TABLE 5-2 


SUMMARY OF LABORATORY CATCHMENT EXPERIMENTS WITH 


UNIFORMLY DISTRIBUTED RAINFALL INPUTS 


Run No. Rainfall Duration Rainfall Intensity Bed Slope 
(sec) (em/hr) So = sin 6 

1 50.0 11.53 0.20791 

2 ils " " 

3 8.8 ‘ 

4 50.0 Velt i 

5 1455 ty - 

6 11.8 . " 

/ 50.0 3.84 . 

8 2500 . " 

9 19.0 " " 
10 50.0 11.47 0.10450 
hi 14.5 u " 
12 TiS uf " 
13 50.0 Ay 42 C 
14 Us) : i 
15 14.5 M " 
16 50.0 SL s 
17 22d " " 
18 18.5 t " 
19 50.0 11.04 0.05117 
20 14.5 i " 
21 Lied : " 
22 50.0 Ton . 
23 Lied 4 4 
24 Vans a " 
25 50.0 3.65 ¢ 
26 30.5 " " 


27 23.9 - d 
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Simulated hydrographs by the state variable model for low rainfall 
intensity and slopes 12° and 2°56', compared to the recorded hydrographs, 
are shown in Figures 5-15 to 5-20. The rest of the simulated events 
(Table 5-2) are shown in Appendix E. 

In most cases, the simulated hydrographs agree very well with the 
recorded events, generally following an "average'' curve indicated by the 
scatter diagram of the recorded hydrographs. The scatter of points 
sometimes becomes larger, or the simulated peak discharge is higher than 
the recorded one, such as can be seen in Figures E-3 and E-19, Appendix 
E. These facts can be partly attributed to errors of recorded data 
which were due to: (1) long sampling interval - about 2 seconds, 

(2) internal noise of the data acquisition system, and (3) manual control 
of the on and off commands for rainfall input. 

The recorded data were smoothed during computer processing and 
then differentiated numerically. Numerical differentiation is sensitive 
to data fluctuation and can become an unstable process. Hence the 
scatter of points. Because of this scatter the comparison between actual 
and simulated hydrographs cannot be carried into fine details. The 
results could be improved by increasing the sampling rate. Unfortunately, 
this was not possible because, as pointed out in Chapter II, the data 
acquisition system malfunctioned during faster rate of sampling. Accuracy 
and error estimation problems are further discussed in Chapter VI. 

The retarding effect of rainfall on sheet flow is particularly 
evident in the simulated hydrograph shown in Figure 5-18. Cessation of 
rainfall caused a sudden increase of discharge. Although this increase 
in discharge could not be confirmed by the experimental data (because 


of their scatter), such phenomena have been observed (Izzard, 1944; 
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Woo and Brater, 1962; Yu and McNown, 1964; Das, 1970; Morgali, 1970; 
Schreiber and Bender, 1972). The increase of discharge does not always 
happen to the same degree. It basically depends on the flow velocity, 
rainfall intensity, and possibly on the surface roughness. Surfaces 
with mild slopes tend to produce greater increase in discharge after 


cessation of rainfall than those with steep slopes. 


5-5-2 Experiments with Temporally Varied Rainfall Input 


Two types of the temporally varied rainfall input were tested: 
(1) ramp function input, and (2) two consecutive pulse function inputs, 
separated by a "no input" time interval. 

Results of the ramp function input experiment are shown in 
Figure 5-21. The simulated hydrograph is a good approximation to the 
laboratory data. 

The two pulse function inputs produced a double peak hydrograph, 


shown in Figure 5-22. The fit of the simulated hydrograph to the 
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experimental data is good. One phenomenon, not observed in the laboratory 


data (again possibly because of the scatter of the data and because of 
too long a sampling interval), is the sudden drop in discharge at the 
beginning of the second pulse input. This is due to increased friction 


factor, caused by the raindrop impact effect. 


5-5-3 Experiments with Spatially Varied Rainfall Input 


The physical phenomenon of surface runoff from a catchment under 
spatially distributed rainfall differs from that under uniform rainfall 
by the difference in distribution of the local inputs and consequently 


by the difference in the time distribution of the output. 
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For experimental purposes two basic patterns of spatially distri- 
buted rainfall input can be distinguished: (1) moving rainstorm with 
uniform rainfall intensity, and (2) stationary rainstorm with non-uniform 
spatial intensity. The former type was studied by Yen and Chow (1968). 
Combination of the two types will result in both spatially and temporally 
varied input. If on such input still another temporal variation is super- 
imposed the result is then the most general case of rainfall input, such 
as most likely occurs in nature. 

From the point of view of runoff simulation, any type of spatially 
varied rainfall input may present a difficult problem. The spatially 
distributed input induces an uneven depth of flow in various parts of 
the test basin, which may eventually lead to formation of waves moving 
over the runoff plane. Simulated outflow from the basin under such 
circumstances, by the state variable model or any lumped parameter 
model, could become significantly distorted. The state variable model 
treats the wave input into a particular reach of the model in the same 
manner as the rainfall input, i.e. as an impulse resulting in an 
immediate response in the form of increased outflow from that reach. 

Such response approximates the actual situation very closely when the 
significant input is in the form of rainfall, but in the case of a wave 
entering the upper part of the reach there is a time lag involved, equal 
to the travel time of the wave through the reach, before there can be 
any apparent response, 

A similar situation arises in the case of receding flow when there 
is no input due to rainfall. The flow is varied gradually and the inputs 


and outputs for a particular reach are of the same order of magnitude. 
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Such gradually varied flow obviously does not present any problem during 
simulation, as can be geen from the very good fit of most of the simulated 
recession hydrographs. 

One can therefore reason that unless there are extreme spatial 
differences in the rainfall intensity, able to induce formation of waves 
with depths several times the magnitude of the flow depth in the rest 
of the basin, the state variable model should provide reasonably good 
approximations to runoff caused by spatially varied rainfall inputs. 

Such extreme conditions, however, are unlikely to occur in overland flow 
situations. 

To examine the effect of spatially distributed rainfall input on 
the simulated hydrograph by the state variable model a stationary 
rainstorm with alternating regions of "rain" and "no rain'' input, as 
shown schematically in Figure 5-23, was generated over the test catch- 
ment. The simulated hydrograph approximates reasonably well the 
experimental data. Unfortunately, the exact shape of the hydrograph, 
especially at the beginning of the rising stage and after cessation of 
rainfall when a sudden increase of discharge occurs in the simulated 
hydrograph, could not be confirmed by the experimental data because of 
the scatter of points. 

Although the detailed form of the response could not be confirmed, 
the experiment demonstrated the ability of the state variable model to 
simulate the hydrograph. It was not attempted to determine to what 
extent, or extremes in spatial distribution of the input, this would hold 
true. It was reasoned that the size of the laboratory catchment ( 2 ft 


by 3 ft) was not sufficient for this purpose. 
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CHAPTER VI 


DISCUSSION OF EXPERIMENTAL RESULTS 


AND SIMULATION METHODS 


6-1 General 

A question that may be raised in connection with the significance 
of the results obtained in the present study relates to the accuracy of 
the method of measuring the outflow rates by numerical differentiation 
of the cumulative outflow. Sufficient accuracy of measured hydrographs 
is necessary to lend significance to the conclusions related to the 
simulating techniques, made through comparison of the measured and 
simulated hydrographs. In this chapter estimates of the accuracy of the 
experimental techniques used in the present study are discussed. 

This chapter also deals with the comparison of the mathematical 
models of surface runoff presented in this study, i.e., the state 
variable model, kinematic wave model, and the instantaneous unit hydro- 
graph model. The aim of the discussion is to investigate the physical 
validity of the basic concepts of each of the above mathematical models. 
Comparison of the output hydrographs of the models simulating the same 
rainfall event serves to depict properties characteristic to each 
mathematical model and to discuss these properties in terms of the 
underlying concepts. 

The last section of this chapter deals with the possible appli- 


cation of the state variable model in simulation of natural watersheds. 
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It is shown that the concept of state and state variables, applied in 
the present study to the elementary case of runoff from an impervious 
surface, is not only useful but perhaps indispensable in the formulation 
of watershed models based on physical principles. This conclusion was 
arrived at by reviewing the present state of the art of digital 


simulation in hydrology. 


6-2 Accuracy of Experimental Techniques 


In the following analysis the dispersion of the computed dis- 
charges is discussed in the case when the number of repetitions in each 
run was seven. The usual procedure is to establish a confidence band 
of the experimental curve given by the average values of the repetitions. 
Because the instantaneous outflow rates were obtained from the measured 
cumulative outflows by numerical differentiation, which is an unstable 
process, the calculation of confidence bands for the instantaneous out- 
flow rate hydrographs must consider possible errors due to differentiation. 
Typical results of calculation of the confidence limits for both 
cumulative and instantaneous outflows are presented in the following 


sections. 


6-2-1 Cumulative Outflow Hydrographs 


Cumulative outflow hydrographs were obtained by recording the 
weight of accumulated water from surface runoff. Accuracy of the 
cumulative outflow hydrographs depends directly on the following two 
factors: (1) the experimental technique, and (2) the overall performance 
of the data recording system, discussed in more detail below. 


The experimental technique of measurement involved a manual 
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operation of a sliding tray by means of which the duration of rainfall 
over the laboratory catchment was regulated. It was estimated that the 
maximum error in timing due to the manual operation was about +0.5 
second. Thus the cumulative outflow hydrographs could be shifted along 
the time axis by a maximum amount of +0.5 second from the true position 
which then results in corresponding differences or errors in values of 
the cumulative outflow. The errors increase with increasing slope of 
the catchment and with increasing rainfall intensity. This is due to 
higher velocities of flow induced by steeper slope and higher intensity, 
which result in larger runoff volumes, and therefore greater errors, 
during the same time interval. Estimated maximum errors in cumulative 
outflows, corresponding to different intensities, are given in Table 6-1. 

The second factor influencing the accuracy of the experimental 
data is a function of the overall performance of the data recording 
system which consisted of the strain gauge transducer that measured 
cumulative weight of water and the data acquisition system. It was 
determined that the error in the reading was +2.0 grams for static weight 
measurement. The error of reading, however, increased to +3.0 grams 
when the weight increase was due to inflowing water from the laboratory 
catchment. The impact of inflowing water into the collecting tank 
caused slight vibrations and waves which were reflected in wider 
fluctuation of the transducer's output signal. 

Estimated maximum errors in the experimental values of cumulative 
outflow are given by the sum of errors due to the two factors discussed 
above. The estimates of maximum errors for three typical rainfall 


intensities are given in Table 6-l. 
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Statistical analyses of the actual test series were also performed. 


Typical results of these computations are presented in Appendix F. 


TABLE 6-1 


ESTIMATED MAXIMUM ERRORS IN CUMULATIVE OUTFLOW 


Rainfall Intensity 
em/hr 115 ieee. 3.9 


Error Due To 
0.5 Sec Time 
Shaft 
(grams) 


I+ 


8.6 260.0 £3 i 


Error Due To 
Recording 
(grams) 


I+ 


oy Af iso gl) 3.0 


Tfotal Error 
(grams) ml HEY) +9.0 +6.0 


In general, at the 95 per cent confidence level, the cumulative outflow 
values determined from seven runs could be considered accurate within 
about 0.9 to 5.5 cc. Figure 6-1(b) shows the 95% confidence bands 

for the series of experiments of run no. 7 (Table 5-2). These results 
are within the maximum values estimated in Table 6-l. 

Figure 6-1(b) also shows the cumulative flow hydrograph predicted 
by the state variable model. It can be seen that the predicted hydrograph 
is within the 95% confidence zone of the experimental hydrograph. The 
close fit of the predicted hydrograph to the experimental values was 
utilized in the following section for estimation of errors due to 


differentiation. 
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6-2-2 Instantaneous Outflow Hydrographs 


Determination of the confidence intervals of the instantaneous 
outflow hydrographs consisted essentially of the following steps: 

(1) Fitting a smooth theoretical curve to the experimental 
cumulative outflow hydrograph. 

(2) Differentiation of the fitted theoretical curve and deter- 
mination of the errors due to numerical differentiation. 

(3) Comparison of the theoretical and experimental differential 
curves and estimation of the confidence interval of the experimental 
curve. 

As pointed out in the previous section, the cumulative outflow 
hydrograph simulated by the state variable model gave a very good fit to 
the experimental curve (Figure 6-1(b)). The simulated hydrograph was 
therefore taken as a theoretical curve fitted to the cumulative outflow. 
The theoretical values of outflow were known at equidistant time intervals 
of 0.1 second. The errors due to differentiation of the theoretical 
curve were determined in the Appendix F, and were found to be less than 
or equal to 0.0007 cc/sec. Since the measured steady outflows were of 
the order of 5 to 18 cc/sec the error in the theoretical curve was 
negligible. 

The theoretical and experimental curves of instantaneous outflow 
rates are plotted in Figure 6-1(a) for comparison. Also plotted is the 
95% confidence interval determined from the experimental data after 
differentiation. Typical results of the confidence interval calculations 
are presented in Appendix F. Since the theoretical curve lies within 
the 95% confidence limits the interval was accepted as valid. The 


curves shown in Figure 6-1 are typical results obtained in this 
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investigation. It can be noted that the experimental hydrograph has 
relatively high values at the beginning of runoff, where they should be 
close to zero. Also the experimental curves tend to be rounded off at 

the time of rain cessation, where there should be a sharp break between 

the end of steady runoff and the beginning of recession. These inaccuracies 
were attributed to the smoothing procedure applied to data prior to 
differentiation. For the same reason it is suspected that the experimental 
values of peak outflow rates were in some cases distorted and the con- 
fidence limits around the peak are unreliable. 

In summary, it was found that with sets of seven consecutive 
experiments the outflow rates could be estimated with 95% confidence 
within approximately 0.1 to 0.4 cc/sec. This represents between 2 to 
8% of the measured steady flow rates. In general the outflow rates for 
various inputs differed from each other by much more than these 
percentages and it was concluded that the method based on the numerical 
differentiation of the cumulative curves was adequate for detecting the 
kinds of responses normally encountered in the investigation. The 
evaluation became uncertain only when the contrasts between responses 
were smaller than about 5%. 

A more efficient means of regulation and control of the rainfall 
input would improve accuracy of reading substantially. Also the high 
sensitivity of the strain gauge transducer was somewhat countered by 


vibrations of the collecting tank. 


6-3 Comparison of the Studied Models 


Simulations of the experimental runs number 13, 14, and 15 


(Table 5-2) were chosen for comparison of the state variable, the 
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kinematic wave, and the instantaneous unit hydrograph models. The 
specific conditions of the runs were: Sp, = 0.1045, rainfall intensity 
Iosu,.42 cn/nr, rainfall duration t =.50.0 sec, 17.5 sec, and 14.5 sec. 
The hydrographs, simulated by the above mathematical models, were plotted 
for each run on a common plot, without the experimental hydrographs. 

The plots show typical results of simulation by the individual models, 
and provide a suitable ground for discussion of the concepts and 


assumptions which form the basis of a particular model. 


6-3-1 The Instantaneous Unit Hydrograph Model 


The basic issue in discussion of the instantaneous unit hydrograph 
model is not the concept itself. The concept of convolution is correct 
and valid, but only for linear processes. Thus, the main concern of the 
investigation was to find whether surface runoff can be approximated by 
the linear unit hydrograph model to such a degree that the method would 
be of practical value. 

It was shown in Chapter IV that two instantaneous unit hydrographs 
of the laboratory catchment, derived from two different runoff events, 
differed considerably in peak values, time to peak, and total time base. 
To facilitate discussion in the present section, another two instantaneous 
unit hydrographs were derived from runs number 13 (i = 7.42 cm/hr) and 
16 (i = 3.91 cm/hr). The derived IUH-1 and IUH-2 are shown in Figures 
6-2 and 6-3 respectively. The large difference between the two IUH's 
indicates a high degree of nonlinearity of the runoff process and 
confirms the results presented in Chapter IV. 

To understand the implications of these results it is convenient 


to recall at this point some of the deductions concerning nonlinear 
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Figure 6-2. Instantaneous unit hydrograph (corresponding to unit 


rainfall intensity i = 7.42 cm/hr) derived from test 13. 
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Figure 6-3. Instantaneous unit hydrograph (corresponding to unit 


rainfall intensity i = 3.91 cm/hr) derived from test 16. 
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Systems, made earlier in this study when dealing with the analysis of a 
nonlinear reservoir. In Chapter III it was shown that a nonlinear pro- 
cess can be approximated by a linearized relationship of the process 
variables around a steady state, and that the approximation is then 
reasonably good only within a narrow range of fluctuations of the process 
variables about the particular steady state. This general conclusion 
can be applied in the case of surface runoff as described below. 

Results of the laboratory catchment experiments indicate that the 
surface runoff is a nonlinear process. The experimental evidence also 
indicates that when deriving the laboratory catchment IUH from measured 
hydrographs generated by rainfalls of different intensities, one obtains 
a number of different IUH. It follows then that there exists an infinite 
number of instantaneous unit hydrographs for a particular catchment, and 
each of these hydrographs is in effect a steady state linear approximation 
(i.e. an approximation valid for a particular value of rainfall intensity 
with infinite duration) of the nonlinear rainfall-runoff relationship. 
Hence, the instantaneous unit hydrograph model is capable of simulating 
only that runoff hydrograph from which the IUH used in the model during 
simulation was derived. When the input amplitude, i.e., the rainfall 
intensity, increasingly deviates from the original "steady state" value, 
the simulation becomes more and more unrealistic. 

Results of the instantaneous unit hydrograph model simulation, 
using the two instantaneous unit hydrographs derived from runs number 
13 and 16 are compared with hydrographs simulated by the kinematic wave 
and state variable models in Figures 6-4, 6-5, and 6-6. Simulated 


hydrographs for the case when rainfall duration was longer than the time 
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to equilibrium (run 13) are shown in Figure 6-4. It should be noted 
that the IUH-1 was derived from the state variable model hydrograph shown 
in the figure because the latter was judged to be a very good fit to the 
experimental data. It is seen that the original state variable model 
hydrograph is oe perfectly by the instantaneous unit hydrograph 
model (using IUH-1) up to the time of rainfall cessation. When the 
IUH-2 was used in the simulation the reproduction was very poor. 

Simulated hydrographs with rainfall durations shorter than the 
time to equilibrium (runs 14 and 15) are shown in Figures 6-5 and 6-6. 
The IUH-1 was used in the simulation. It is seen that the original 
(state variable) hydrographs are reproduced again correctly up to the 
time of rainfall cessation. Simulation of the recession curves is 
unsuccessful in all cases. The cause of this failure can be found in 
the shape of the IUH, which is not suited for reproduction of the 
recession curves by means of a convolution integral, and the reason can 
be attributed to the nonlinearity of the runoff process. 

In summary, the performance of the instantaneous unit hydrograph 
model was found unsatisfactory. The main reason for the poor simulation 


results was found to be the incorrect theoretical assumption of linearity. 


6-3-2 Kinematic Wave and State Variable Models 

Results of the present study have shown that the kinematic wave 
and state variable models were able to simulate reasonably well rainfall- 
runoff events within the laboratory conditions encountered in the 
investigation. In the following the properties of the kinematic wave 
model by the way of comparison with the state variable model are discussed. 


This method of discussion was selected primarily because the principle of 
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the state variable model is simple to the degree that many flow 
situations can be easily visualized. 

The kinematic wave model as well as the state variable model of 
surface runoff are based on the assumption that the discharge of over- 
land flow can be described with a reasonable accuracy by a uniform flow 
formula. The two models, however, differ in the mathematical manipulation 
of the equation of continuity of the unsteady spatially varied flow. 

In the mathematical development of the kinematic wave model the 
equation of continuity was written as a total derivative of depth of 
flow [Equations (4-13), (4-16), and (4-18)]. Validity of the new 
equation thus derived was based on the assumption that discharge is a 
function of depth only, and was restricted to a new set of co-ordinates. 
Analytical solution of the surface runoff problem then followed from 
the application of the method of characteristics. 

The state variable model is based on the assumption that the 
partial differential equation of continuity can be approximated by a 
finite number of simultaneous ordinary differential equations of con- 
tinuity. The set of such equations can be written in the form of state 
equations. It was shown in Chapter V that as the number of equations in 
the set increases to infinity, the set approaches the partial differential 
equation of continuity of the flow in the limit. This property makes the 
state variable model physically meaningful, and suitable to represent 
the actual flow very closely. As a result, the state variable model was 
used in the following to estimate the validity of some of the relations 
given by the kinematic wave model. 


Two characteristic properties of the kinematic wave model discussed 
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in this section are the following: (1) the absence of a point of 
inflection on the rising limb of a hydrograph, and (2) the steady rate 
of receding piterde for the period of time tp, when to<ts. 

The first of the above properties is illustrated in Figure 6-4. 
The absence of a point of inflection is the result of the kinematic wave 
relationship specifying the depth of flow during the rising stage as 
y = it (or dy/dt = i) given by Equation (4-19). The peak of the 
simulated hydrograph is then always sharp crested as illustrated in 
Figures 6-4 and 6-6. Such a shape of hydrographs has never been observed 
in actual outflows and the validity of the relationship y = it is 
therefore questionable. 

In the following, the elementary case of runoff, i.e., runoff 
from a uniformly sloping plane under uniformly distributed rainfall, is 
analysed by means of the state variable approach to assess approximately 
the validity of the relation given by y = it for tx tg. 

Consider for a moment an infinite runoff plane. The state 
variable model simulating runoff from this plane would consist of an 
infinite number of identical elements in series, basically operating as 
a series of interacting reservoirs (Chapter V). Since the rainfall is 
uniformly distributed, each of the elements receives exactly the same 
rainfall input, and the depth in all elements rises {iis Hah we the 
first time increment) to the same depth. Consequently, the outflows 
from the elements are identical. Therefore, during the following time 
increments the surface outflow and inflow components of the continuity 
equations of the elements are all equal, and the storage components are 


made up only of the precipitation input, i.e. y = it. 
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Consider now a runoff plane of a finite length. The first 
(upstream) element in the series does not receive any surface inflow 
but only the precipitation input. Because the outflow from the element, 
q, is a function of the third power of the depth, y, i.e., q < y3, the 
response of the element to a step function input has the form of an S- 
curve. Let Q) be the rate of outflow from the first element corresponding 
to the point of inflection of the S-shaped outflow curve from the element. 
The S-shaped outflow function is routed through the rest of the series, 
which disturbs the equality of inflow and outflow, described for the 
case of infinite series. Consequently, the process results in differentiation 
of depths in the individual elements of the series. The outflow hydrograph 
simulated by the kinematic wave model corresponds to the outflow from 
the last element in the series representing the state variable model. 
It is reasoned that for the last element in the series the relationship 
given by y = it remains valid for the time required to route the discharge 
Q, through the entire series. After the discharge Q)] reached the n-th 
and last element of the series, the rate of increase of the surface 
inflow into the n-th element begins to decrease, thus marking the point 
of inflection of the total outflow hydrograph. The total outflow then 
approaches assymptotically a steady state rate of outflow. Table 6-2 
illustrates the described process on an example of the run number 4 
(Table 5-2), simulated by the state variable model with n = 7. The rows 
of the table give the simulated water surface profiles at small increments 
of time. 

The discussion presented above leads to the conclusion that both 


the kinematic wave and state variable models, simulating a particular 
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runoff event, will compute the same values of discharge up to the time 

of occurrence of a point of inflection on the hydrograph simulated by 

the state variable model. For illustration, hydrographs of run number 13 
(Table 5-2), simulated by the kinematic wave and the state variable 
models are plotted in Figure 6-4. It can be seen that the hydrographs 

do not differ in their shape significantly up to the time of occurrence 
of a point of inflection on the hydrograph simulated by the state 
variable model. The faster rate of rising of the kinematic wave hydro- 
graph is due to neglected raindrop impact effect. Only after the time 

to the inflex point do the two hydrographs differ in their shape signifi- 
cantly. 

The second property of the kinematic wave model listed previously 
is typical for the case when the duration of rainfall, to, is shorter 
than the time to equilibrium, tg. In this case, after cessation of 
rainfall the discharge remains constant for the time period tp. A 
typical outflow hydrograph is shown in Figure 6-6. Considering this 
case from the point of the state variable model, the constant outflow 
could occur only if the depth of flow, at the time of rain cessation, 
is the same in two or more elements of the series. Referring to the 
previous analysis of flow by the state variable approach, this means 
that in order for the receding flow to remain constant for a certain 
period of time, tp, the rain must cease before the occurrence of a point 
of inflection on the hydrograph of total outflow, caused by rain of 
duration to > tg. The flow during the time period tp can be visualized 
as propagation of a negative wave through the remaining elements of the 


series. Clearly, the time period tp depends on the chosen number n of 
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the elements in the series, and on the time interval used in the calcu- 
lation. These two factors must be determined through experience, and 
will require the use of an empirical coefficient in the outflow equation 
(Equation (5=51)). Such a coefficient should be determined by a parameter 
optimization technique. In actual flow the time period tp» if it exists 
at all, most likely depends on the slope of the runoff plane and on the 
energy losses due to friction. The faster the flow the shorter the tp 
will be. Hydrographs simulated by the state variable model (Chapter IV 
and Appendix E) do not show any time period tp of a significant duration 
throughout which the flow would remain strictly constant. However, some 
of the recession curves are S-shaped, and on very mild slopes the rate 
of decrease of discharge is so small for an appreciable duration that it 
is almost constant. A typical case is shown in Appendix E, Figures E-15 
and E-21. 

Recession curves of equilibrium hydrographs, i.e., the case when 
to > ts, simulated by the kinematic wave and state variable models are 
practically identical, as illustrated in Figure 6-4. The state variable 
model recession curve is only shifted horizontally due to increased 


runoff caused by decreased friction factor at the time of rain cessation. 


6-4 Application of the State Variable Approach in Hydrologic Simulation 


With the rapidly increasing rate of urban development in the 
present technological society there is an ever growing need for better 
and more effective methods of water resources management, flood regulation, 
streamflow forecasting, pollution control, and power regulation. All 
these needs require development of more reliable techniques of hydrologic 


forecasting. Taking as an example the problem of flood forecasting, the 
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simplest method of obtaining a flood forecast at a downstream point would 
consist of monitoring river and reservoir levels upstream and routing 
the flood hydrograph downstream. In many instances, however, the travel 
time of the flood waves through the channel reaches “ig short. Thus, 
there would be insufficient time to prepare and disseminate the forecast 
before the onset of flooding. Any effective forecasting system must 
give a much earlier warning and this implies the necessity of basing 
forecasts on measured or forecasted rainfall which must then be routed 
mathematically through the watershed as well as the river channel and 
reservoir system (City of Calgary Flood Study, 1972). 

The transformation of measured or forecasted hydrometeorological 
data into a flood forecast may be achieved by means of a watershed 


"models'' the behaviour of the watershed. A 


simulation system which 
brief review of the state of art of hydrologic simulation is given in 


the following. 


6-4-1 Traditional Approach 


In the simplest case of a fairly small uniform watershed in its 
natural state affected by fairly uniform rain only, one of the traditional 
and most common flood forecasting procedures would involve, for example: 

(1) Determination of the storm rainfall on the watershed and its 
duration by plotting all measured rainfall data on a map, sketching in 
the isohyets and planimetering them to obtain the estimated watershed 
rainfall. 

(2) Determination of the current antecedent precipitation index 
(API). 


(3) Determination of the storm runoff volume from a coaxial 
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rainfall-runoff diagram using watershed rainfall, its duration, API, and 
time of year. 

(4) Determination of the flood hydrograph at a gauged point for 
which a unit hydrograph is available. 

(5) Determination of the flood hydrograph at the forecast point 
by routing the flood hydrograph determined through the intervening reach 
by the Muskingum or other routing procedures. 

The traditional approach to flood forecasting is deceptively 
simple in appearance but involves considerable effort in practice when 
complications arise from the non-uniformity of meteorological and surface 
conditions and the resulting hydrologic complexity. Aside from these 
considerations, the method involves principles which are known to be 
incorrect. The coaxial diagram is a very crude model of the runoff 
process and involves statistically fitted curves of best fit to data 
which are widely scattered. The principle of the unit hydrograph has been 
questioned for many years. The investigation in the present study has 
revealed that the unit hydrograph concept is theoretically unsound. 
Most of the traditional methods therefore, are not suitable for operational 


flood forecasting. 


6-4-2 Modern Approach - Digital Hydrologic Simulation 


The advent of high speed digital computers has removed the 
traditional limitation of calculating speed. The object of today's 
research is to develop a general system of quantitative analysis for 
hydrologic regimes by establishing continuous mathematical relationships 
between elements of the hydrologic cycle. The operation of these 


mathematical relationships is observed and improved by using digital 
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computers to carry the calculation forward in time. Precipitation and 
potential evapotranspiration are the basic inputs for these digital 
models and actual evapotranspiration, streamflow, and soil moisture 
levels are generally obtained as outputs. Calculations are made on 
selected time intervals and are carried continuously, whether or not 
precipitation is occurring, to simulate the entire spectrum of watershed 
behaviour. 

A review of the literature has indicated that there are a dozen 
watershed simulation systems presently in existence. Most of then, 
however, are either in the development stage or are specifically intended 
for some purpose. An example of such "dedicated" model is the streamflow 
simulation and Reservoir Routing Model developed by the U.S. Corps of 
Engineers and the Portland River Forecast Centre for use on the Columbia 
River Basin. It is essentially a flood routing model which lacks adequate 
simulation of the land surface runoff process. 

Perhaps the best known digital hydrologic model is the Stanford 
Watershed Model (Crawford, Linsley, 1966). The authors have developed 
the model over the last 12 years into a powerful program known as the 
Hydrologic Simulation Program (HSP). HSP programming is divided into 
three load modules designated Library, Lands, and Channels. The Library 
module handles the hydro-meteorological data files. The Lands module 
simulates land surface processes such as infiltration, evaporation, soil 
moisture storage and overland flow. The Channels module simulates depth 
and velocity of flows in the stream channel system, and spillway dis- 
charges and reservoir levels. 

Two main steps are required to bring HSP into operation. These 


are "installation" and "calibration". Installation comprises establishing 
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the basic mode of operation, selecting the data network stations which 
will provide input to the program, subdividing the watershed into 
"segments" and the channel network into "sections", and establishing rule 
curves and operating policies for each reservoir. 

Calibration refers to the selection of the parameter values to be 
used in the program and their successive adjustment during trial runs to 
the point where the simulated record fits the observed records as closely 
as possible. A minimum of five years data is required for this purpose. 
This program is extremely versatile and practical, and finds applications 


in flood forecasting and water resources management. 


6-4-3 Potential of the State Variable Approach 


The state variable approach to system analysis appears to be 
suited to digital hydrologic simulation. The advantages of the method 
and its practical application are discussed in the following. 

Because of the almost unlimited diversity of the boundary and 
initial conditions characteristic to the hydrological processes occurring 
in natural watersheds, pure distributed hydrologic models are not 
practical. Instead, the behaviour of a watershed must be simulated by 
lumped parameter models. The state variable method is in many instances 
the most efficient method of lumped parameter systems analysis. The 
method is extremely versatile and enables the mathematical characterization 
of lumped systems in a compact matrix form suitable for digital computer 
calculations. The concept of state itself is invaluable for development 
of a hydrologic model that operates using physically relevant components. 

Application of the state variable method in simulating surface 


runoff under laboratory conditions was investigated in the present study. 
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To apply the method to natural watershed conditions, the state equations 
must be expanded to include evapotranspiration, infiltration and surface 
detention terms. Equations relating these terms in the are and output 
equations must be carefully considered, and will involve empirical 
coefficients determined by parameter optimization. 

Continuous simulation of inflow into the channel network of a 
watershed requires simulation of the subsurface components of this 
inflow. The state variable model of the subsurface flow would not differ 
in its structure from the surface runoff model. The infiltration amounts, 
which are outputs of the surface flow model, will become inputs into the 
subsurface flow model. It is anticipated that the density of sub- 
division, or the number of elements required for subsurface flow 
simulation, would be less than for the surface flow model because of 
smaller flow velocities within the porous media. A block diagram of the 
state variable model simulating the lateral inflow into the watershed 
channel network is shown in Figure 6-7. The subsurface flow model can 
be further subdivided into two subsystems simulating the interflow and 
groundwater flow respectively. The channel flow can be simulated in 
principle by the state variable model, however, its application to 
large velocities and depths of flow in the unsteady channel flow would 
require further investigation. Various other methods of channel flow 
routing, however, which have been used successfully are available, such 
as the kinematic wave or Muskingum methods, 

In the practical application of the state variable model the water- 
shed will be divided, for the purpose of installation and calibration and 


subsequent operation of the model, into a number of segments and the 
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channel network of the watershed into a number of reaches. A definition 
sketch of the watershed components used in the model is shown in Figure 
6-8. Each segment and each reach has its characteristic parameter values 
used in the state and output equations of the model. For example, the 
state variable model of surface runoff developed in the present study 
contains only one parameter, a, [Equations (5-60) and (5-61)], which 
depends on the discharge formula and the time and space increments used 
in the calculation. State variable model simulating both surface and 
subsurface flow will contain more parameters to account Peer embaureace 
flow conditions, infiltration, evaporation, etc. The parameters will be 
estimated from existing information on the geology, topography, vegetation, 
channel cross-sections, etc. Almost all this information can be derived 
from existing maps, aerial photographs, bridge surveys, dam surveys, etc. 
The estimated parameters will be optimized to provide best fit to 
observed hydrographs. 

Surface runoff from each segment will be calculated by the state 
variable model of surface flow fitted to each segment individually. As 
mentioned previously, the subsurface flow model of each segment will 
likely consist of a smaller number of elements than the corresponding 
surface flow model, i.e. each element of the subsurface flow model will 
receive inputs in the form of infiltration from several elements of the 
surface flow model. 

The general state variable model will continuously calculate the 
streamflow from rainfall in each successive watershed segment, and from 
flows measured or calculated at upstream flowpoints. All calculations 


will be carried independently for each segment so that areal variations 
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IN A CATCHMENT SEGMENT BY 
STATE VARIABLE MODEL 


tte Gaier 2S 6 ams @ a a 


SEGMENT OF THE 
CAT CHMENT 


CATCHMENT BOUNDARY 


Figure 6-8. Definition sketch for application of the state variable 
model in simulating surface runoff from a natural watershed. 
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in rainfall and topographic features can be represented. There will be 
no limit to the number of segments that can be used. 

In summary, the concept of state and the state variable method 
form a skeleton of a practical and physically relevant hydrologic model. 
The state variable model can reasonably well represent spatial and 
temporal variations of rainfall and topography, and is suitable for 
continuous simulation by a digital computer. The structure of the model 
allows for improvement of the performance of the model by: (1) changing 
the inter-relationships between individual components, such as the 
infiltration, evapotranspiration, etc., as they become better understood, 
(2) optimizing parameters during the simulation, and (3) increasing the 
number of segments or elements into which the watershed has been sub- 
divided. The property (3) gives the state variable model considerable 
flexibility since it can be easily applied to different watersheds. For 
the same reason the state variable model has also potential for use in 


simulation of ungauged watersheds. 
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CHAPTER VII 
CONCLUSIONS AND RECOMMENDATIONS 


7-1 Conclusions 

The study shows that unsteady spatially varied overland flow 
from an impervious plane can be simulated by a number of interacting 
elements operating as reaches with uniform flow subjected to impulse 
inputs. Mathematically the simulation model is conveniently described 
by the state space method of system representation. 

The concept of state enables expression of the antecedent con- 
ditions, i.e., to determine the state of a watershed, by means of a set 
of state variables which have a clear physical interpretation, such as 
the depth of flow. The mathematical form of the model is suitable for 
computer solution. 

The state variable model developed in the present study contains 
only one parameter, a, which depends primarily on the discharge formula 
used (for laminar or turbulent flow) and halomas the time and space 
increments used in the calculation. No attempt has been made to derive 
a criterion for the choice of the time and space increments because of 
the experimental limitations imposed by the size of the laboratory 
catchment. It was found, however, that the calculated hydrographs did 
not differ significantly for time increments of 0.1 sec and 0.01 sec 
respectively. With respect to the number of segments used in simulation 
it was found that increasing the number of segments resulted in a 


faster response to a given input. Theoretically, there could be an 
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infinite number of segments. For practical reasons, however, it was 
concluded that some judgement must be used to find an optimum number of 
segments which would be consistent chars accuracy of other input 
parameters, computer storage availability, etc. The use of a finite 
number of segments may require empirical evaluation of a correction 
factor applied to the discharge (output) equation. 

The study shows that rainfall affects the mechanics of overland 
flow on a hydraulically smooth boundary. The rainfall intensity is the 
strongest single parameter affecting the flow. The relationship f;=Cr/Re 
describes adequately the flow resistance of overland flow under rainfall 
within the range of experimental conditions covered in the investigation, 
i.e., Re <30, U < 10 cm/sec, 3 < i< 12 cm/hr, and for the galvanized 
wrought iron runoff surface treated with hydrochloric acid, having 
Manning's n approximately equal to 0.015. The coefficient Cr is a function 
of rainfall intensity, Reynold's number and the velocity of flow. 

Simulated hydrographs by the state variable model showed increased 
outflow rates after rain cessation in the form of a sudden peak, This 
increase in discharge was attributed to the change in the friction factor 
fr; this coefficient is a function of rainfall intensity. 

The state variable model is capable of defining overland flow for 
both cases of uniformly distributed rainfall and temporally and spatially 
varied rainfall. The results showed that the calculated hydrographs 
were in good agreement with the experimental data. The fit of the 
simulated hydrographs to the experimental data was slightly better for 
milder slopes and lower rainfall intensities. 

The experimental results showed the surface runoff from an 


impervious surface to be highly nonlinear. The apparent time to 
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equilibrium, a measure of the time required to reach steady state out- 
flow, for fixed physical parameters of the catchment (length and roughness) 
was found to be a function of rainfall intensity, contrary to the unit 
hydrograph theory which assumes a constant time to equilibrium. 

Instantaneous unit hydrographs (IUH) derived by differentiation 
from observed experimental equilibrium hydrographs (S-curves) differed 
from each other by as much as 50 per cent in peak values, time to peak, 
and total time base. Accordingly, computed hydrographs by the 
instantaneous unit hydrograph model were unrealistic. It was concluded 
that there exists an infinite number of IUH for for a given catchment, 
and each represents a linear approximation of the rainfall-runoff process 
for a given rainfall intensity. As the rainfall intensity changes the 
particular linear approximation becomes progressively worse. 

The kinematic wave model of surface runoff was found to predict 
the experimental hydrographs quite satisfactory. The model was able to 
respond correctly to changes in rainfall intensity and slope of the run- 
off plane. However, the model was not suitable to account for temporal 
and spatial variations of rainfall. It also did not include a factor 
for the raindrop impact effect. Hydrographs simulated by the kinematic 
wave model are sharp crested owing to the lack of a point of inflection 
on the rising limb of the hydrograph. It was concluded that because of 
this peculiarity the time to equilibrium, predicted by the kinematic 


wave model would be generally shorter than the observed time. 


7-2 Recommendations 


From the results of the foregoing study the following recommendations 


can be made: 
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(1) Surface runoff experiments should be carried out on a larger 
laboratory catchment. A procedure employing direct outflow rate measure- 
ment and a more efficient means of regulation and control of the input 
could reduce the number of consecutive runs necessary and also produce 
a better definition of details. 

(2) Detailed experimental studies are needed to examine the role 
of the rainfall impact effect and surface roughness in overland flow. 

The transition range of overland flow under rainfall from laminar to 
turbulent need to be investigated in order that a proper discharge 
formula can be selected. 

(3) Experiments should be carried out to simulate surface runoff 
from urban areas by the state variable model and to determine criteria 
for finding an optimum number of segments into which a simulated catch- 
ment need be subdivided. 

(4) The results describing the scope of existing experimental 
conditions and performance of the state variable model should be carefully 
considered in the design of a future general hydrologic model. Application 
of the concepts of state and state variables in general hydrologic 
simulation will result in a worthwhile expansion of the theory and 


experimental data presented herein. 
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STRAIN — RESISTANCE RELATIONSHIP 
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A-1 Load Beam: Stress — Strain Relationship 


Consider a beam with two symmetrically acting forces of equal 
magnitude P, as shown in Figure A-l. It is seen from the stress diagram 


shown in Figure A-2 that 


Oy — gamax _ 

Sy Me Baa k (A-1) 
or 

Oy = ky (A-2) 


where oy is the normal stress at the distance y from the neutral axis. 
With respect to the cross-section diagram shown in Figure A-2 it 
is seen that the resultant normal force is 


Rn = { Gy dA = k i y da =0 (A-3) 
area area 


and the resultant internal moment M is 


m= [ yoyat=k J y2aa-ut (A-4) 
area area 


where I denotes the second moment of the cross-sectional area about the 


horizontal centroidal axis. Thus, 


(A-5) 


HS 


The stress - strain relationship for a uniaxial state of stress 
is expressed by Hooke's Law 

o = Ee (A-6) 
where E is the modulus of elasticity and ¢« is the strain. Combining 


Equations (A-2) and (A-6) gives 


e = Et (A-7) 
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Figure A-l1. Definition sketch of a beam in bending. 


cross - section normal stres diagram 


Figure A-2. Stress distribution diagram. 
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A4 
= Ply s 
eek | (A-8) 
and 
€ i« Pp | (A-9) 


It is thus shown that in the case of elastic bending of a beam with 
symmetrical sections the strain is linearly proportional to the applied 


force P. 


A-2 Strain Gauge: Strain - Resistance Relationship 


The mechanical strain e« in a beam will give rise to consequent 
electrical resistance change, AR, in a strain gauge mounted on the beam 
in accordance with the relationship 


aoe 
R 


where K is a gauge factor. Gauge factor (G.F.) is a property of the 


eK (A-10) 


strain gauge which relates output resistance change to input mechanical 


strain and is defined as 


C.F unit resistance change in the strain gauge 


~ mechanical strain in the specimen beneath the gauge 


A-3 Sensitivity of Load Transducer 


When the instrument with which the strain gauge circuit output is 
measured is a voltage-sensing device, the circuit output sensitivity may 


be expressed according to Stein (1962) as 


ye deen (A-11) 
Also 

Vel Rok (A-12) 
and 

V=iR (A-13) 
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Thus, : 

Sean Rime oy (A-14) 
where | 

Sy = system voltage sensitivity in units of microvolts circuit 

output per microstrain in the speciment beneath the gauge 

V = voltage output from the transducer at reference applied P 

v = incremental voltage output produced by the action of P (or e€) 

i = current in amperes passing through the gauge 


Rg = gauge resistance in Ohms 

Nc = efficiency of the circuit of which the gauge is part 

Nm = efficiency of signal transfer from the circuit to the voltage 

measuring instrument 

It can be concluded, by examining Equation (A-14), that the 
sensitivity of a voltage-sensing device used for strain gauge transducers 
can be increased, for a given applied current i, by using strain gauges 
with higher resistance Rg and higher gauge factor K. 


For further reference see Byars (1969) and Stein (1962). 
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APPENDIX B 
DERIVATION OF A UNIT IMPULSE RESPONSE OF 


(1) SINGLE RESERVOIR, (2) SERIES OF RESERVOIRS 
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B-1l Unit Impulse Response of a Single Reservoir 


The system under consideration is shown in Figure B-l. The 
relationship between storage, S, and discharge, q, is linear and has the 


form 


q"_ (B=) 


where k is a constant. Equation of continuity of the reservoir is 


as. 

Ae <Bt57 4 (B-2) 
or 

daae” Gs 

1 ge a (B-3) 


where i is the volumetric rate of inflow. 


S 
q = apa 
k = constant 
= storage 


Donne * AA 
RA AR RA AA AAA AAAS AA | 


Figure B-1. Definition sketch of a single linear reservoir. 
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When the system described by Equation (B-3) operates at a steady 


state, Equation (B-3) becomes 


ase 
lie + saad 
and 
S 
vee —- =40 (Bes) 


where the asterisk denotes steady state values. Subtracting Equation 


(B-4) and (B-5) from Equation (B-3) gives 


d(S=S . 

a a CEI) -< (S-Sy) (B-6) 
or 

dZ . Z x) 

dipencnkek nee 


where Z = S-Sx and I = i-ix are deviation variables. 
Taking Laplace transform of Equation (B-7) gives a transfer 


function between input and storage 


ZAR ae b 
Ti) Sa (B-8) 
k 


The transfer function relating outflow and inflow can be obtained 


as follows. The discharge-storage relationship at a steady state is 


which subtracted from 
q= . (B-10) - 
gives 


1 
q-qx = , (S-Sx) 


Substitution of the deviation variables, Q = q-qx and Z = S-S, gives 


(B-11) 
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Laplace transform of Equation (B-11) is 


S08 Radel 
Te (B-12) 


Multiplying Equation (B-8) by Equation (B-12) gives the desired transfer 


function between input and output 


Q(s) I 
rye FT sek 
For a unit impulse input I(,) = 1 the transfer function given by 


Equation (B-13) becomes 


1 
Q(s) * ks + l (B-14) 
Inversion of Equation (B-14) into the time domain gives the unit 
impulse response of a single reservoir 


u(t) = Q(t) =< e-t/k | (B-15) 


Response of the system to an input i(t) can be found by application 


of the convolution integral 
iz 
q(t) = f i(t) u(t-t)dt (B-16) 
O 


For an input of a constant magnitude i one obtains 


t pines 
q(t) =i J de k dr (B-17) 
eke 
\ 
which integrates into 
a(t) = ife-t/k et/k]™ (B-18) 
or 
q(t) = i(1 - e7t/ky (B+19) 
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B-2 Unit Impulse Response of a Series of Reservoirs 

The system under consideration is shown in Figure B-2. Transfer 
function of a system consisting of n identical linear non-interacting 
reservoirs can be derived by induction. The transfer function between 
input and output for each tank is given by Equation (B-13) and the input 
into the second, third, etc., reservoir is the output from the first, 
second, etc., reservoir. Thus, the transfer function relating the output 


trom the nth reseryoir,7a,, to, the input 1 is 


6 Qnls) © 1 , = 
a eee See 


Figure B-2, Definition sketch of a series of linear non-interacting 


reservoirs. 
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The unit impulse response of the system can be derived conveniently by 
induction. Consider a system consisting of two reservoirs in series. 
The transfer function of the system is 
Q2(s) 1 2 
CU eer Pe ( otcesclal B-21 
I(s) oe as 7% ( ) 
Substituting I(s) = 1 for a unit impulse gives 
1 1 
ee ay ar oe 
2 ke) (oe =)? , (B-22) 
Applying partial fraction technique yields 
1 A B 
Q2(s) = Dipselinecy terpecty (B-23) 
aes oe 
(erie) (sit D 
Inversion into time domain gives 
Qo(t) = e-t/k(At + B) (B-24) 
The constants A and B can be determined from combined Equations (B-22) 
and (B-23) 
a Se (B+25) 


i 
k 
Multiplying Equation (B-25) by a factor (s + 2)? gives A = 1, and 


= tS 
(s + 74 Gr s + 


differentiation of Equation (B-25) with respect to s gives B = 0. Thus, 


the solution of Equation (B-24) is 


Qo(t) = 3 e-t/k (B~26) 


Similarly, for a system consisting of three reservoirs in series, 
one obtains 
iL le 
03 (3) 850 ee (Ba2 7) 
; ke (s + 2)3 
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B7 
which can be expanded into: 
1 A B C 
Q3(s) = =o ee + Perr le + Oat. (B-28) 
+2 += += 
(s > (s 7 st 
and inversion into time-domain gives 
ag(t) = LL e-t/kfat2 + Be + c] (B-29) 
eo 2! 
Rewriting Equations (B-27) and (B-28) as one equation gives 
1 A B 
SALE aa bas sat a a (B-30) 
Multiplying Equation (B-30) by the factor (s + 2)3 gives A= 1, and 
differentiating twice with respect to s gives B = 0 and C = 0. Thus, 
the solution for three reservoirs in series is 
Gy = ale | (B-31) 
Q3 wan 24 = 


By induction the formula for a unit impulse response of a system 


consisting of n linear non-interacting tanks is derived as 


— acl if Sev x 
Qn(t) = a Geel ent / (B-32) 
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APPENDIX C 


RAINDROP IMPACT EFFECT ANALYSES - 


TABLE OF MEASURED AND CALCULATED 


FLOW PARAMETERS 


198 


sy 


eq et SK, 6 j ae 7 bi 
i ’ : vi Pa ae 

Teta f oe bound vy bay Pe Pe at , ’ 
| oa 3 ms 


" i coe il ie one 


es ed “waa 


rien 
a oy 


“i 
+1 


fil +. 
ae 


i nh e i Riek mY 


tok iy 
j . ay - 
1 i ‘ iy 7 : 
; } ; H iain ' ah F 
<a ; ni io Th 
> a ae 
' - 
saa ; oy, ane ees 
} i tf 
aes nH 


RY th Re ee) 
| ne | ati | 


ema Wry" ia 


I 
T/RE*U 
Q 

RE 

lf 


U 


LIST OF SYMBGLS 
= CR=-24 EQe(5~-88) 
= ALPHA*B*(H**3)/Q EQe(5-86) 
= FRICTION FACTOR OF LAMINAR FLOW UNDER 

RAINFALL» FF=CR/RE EQ@e(5S-76) 


FROUDE NUMBER 
EQUIVALENT DEPTH CF UNIFORM FLOW, 
H=ST/(B¥X)« WHERE B AND X ARE THE 
WIOTH AND LCERTH OF THE RUNOFF PLANE 
RAINFALL INTENSITY 
CIMENSIONLESS PARAMETER EQe(5-90) 
STEADY OISCHARGE FROM THE RUNOFF PLANE (CC/SEC) 
REYNOLO NUMBER 
AMOUNT OF DETAINED WATER ON THE RUNOFF PLANE 
DURING STEADY CISCHARGE Q (CC) 
MEAN FLOW VELOCITY IN THE X=DIRECTION (CM/SEC) 
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Program Name: MAIN (including subroutines STATE and WAVE) 


Description: MAIN is designed to perform the following operations: 

(1) Derivation of IUH (MAIN) 

(2) Simulation of surface runoff hydrographs by: 

(a) IUH model (MAIN) 
(b) State variable model (subroutine STATE) 
(c) Kinematic wave model (subroutine WAVE) 

Output from MAIN includes printed tables of time and discharge 
coordinates (1 sec time increment), and Cal Comp plots of derived IUH 
and simulated hydrographs. Subroutine CGPL (available at the Computing 
Center, University of Alberta) was used for plotting. 

MAIN generates an S-hydrograph (using state variable model) of 
surface runoff from a plane of given slope and size, and for the given 
rainfall duration and intensity (alternatively time and discharge 
coordinates of S-hydrograph could be supplied). IUH is derived from 
smoothed S-hydrograph (subroutine SE15) by numerical differentiation 
(subroutine DET3). Subsequently the S-hydrograph (or any other hydro- 
graph) is reproduced (simulated) by convolution using the derived IUH. 
Subroutines SE15 and DET3 are available from standard IBM subroutine 
package. 

Subroutines STATE and WAVE simulate surface runoff hydrographs 
for the same conditions (duration and intensity of rainfall, slope and 


size of runoff plane) as specified in MAIN program. 
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LIST OF SYMBOLS 


A = SPACE INCREMENT 

AA = AOMISSABLE ERROR OF ITERATION 

ALP = PARAMETER ALPHA 

B /= WIDTH OF THE RUNOFF PLANE 

CH = DURATION OF SIMULATED RAINFALL 

C(I) = FRICTION FACTOR 

EC = UPPER LIMIT OF C(I) 

HCI) = DEPTH OF FLOW IN THE I-TH ELEMENT 
HO = DEPTH OF FLOW AT B 

K = TOTAL TIME OF SIMULATION (10*SEC) 

L = LENGTH OF THE RUNOFF PLANE 

M = ORDER OF THE STATE VARIABLE MODEL 
ND = TOTAL TIME OF THE IUH BASE (10*SEC) 
P = RAINFALL INTENSITY (CM/SEC) 

PPP = INTENSITY OF SIMULATED RAINFALL 

Q = DISCHARGE FROM THE RUNOFF PLANE (CC/SEC) 


Q(J) = ORDINATE OF SIMULATED HYDROGRAPH (CC/SEC) 
QUCi)= OROINATE OF THE IUH (CC/SEC) 
RCI) = DISCHARGE FROM THE I-TH ELEMENT 


= = SLOPE OF THE RUNOFF PLANE 

T = TIME (SEC) 

T = TIME INCREMENT (SEC) 

T,0 = DURATION OF RAINFALL (SEC) 

TOT = TOTAL TIME OF SIMULATION (SEC) 

TS = TIME TO STEADY STATE RUNOFF (SEC) 
TT, = TIME INCREMENT (SEC) 


TT(J)= ABSCISSA OF SIMULATED HYDROGRAPH (SEC) 


MAIN PROGRAM 
IUH GENERATION AND CONVOLUTION 


REAL A(10)2B(10).H(10).B8E( 25) «RR 25) »R(10) »Q(800)5TI(400) >» 
IDERY (10) +DH(10) »sALPHA( 20) » Y( 40) 52040) »U( 400) » QU(100) 2 TUC100) » 
1P0OP(200s9) »Z2Z(400)4C(30)3QS(200) s TA( 800) .UQ( 400) 

CCMMON MKsSeALPHAsCHsePPPelLKe IVA 

READ(S5s33) (ALPHACI)».I=1,.20) 

33 FORMAT(20A4) 

IVA=1 

MK=129 
DO 350 JIM=1+s2 

M=7 

G=981. 

V=.0093 

S=0-1045 

TT=e1 

00 60 I=15M 

ACT )=91264*61-/M 

60 CONTINUE 
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ALP=8.*G*S/V 
FC=3le 
DO 2 J=1>M 
C(JI=FC 
H(J)=0- 
2 CONTINUE 
P=0.0012 
IFC JIMeEQe2) P=02-002060838 
PRE=P 
N=1 
It=1 
T=0. 
K=400 
DO 200 L=1>sK 
DO 50 I=15M 
R(I)=ALP*61-¢* (HC I)**32.)/C( 1) 
IFCRC1I)-EQe0-) GO TO 50 
QC=R(1I)7616 
ABC=2 eS¥*¥ALOG( (1000006 *P*¥V*XECI))7(01264*QC*QC) ) 
C(I) =24-.+ABC 
50 CONT INVE 
DO S51 f[=1.5M 
IF(PeEQeO-) C(I)=24. 
RCI V=ALP*61 -*¥ CHC I) *¥*32)/CC1) 
51 CONTINUE 
Q(L)=RCM) 
TACL)=T 
DERY (1) =P-R(1)/7A(1) 
DO 70 I=23M 
DERY CI) =P4CRCI-1)-RCIDISACT) 
70 CONTINUE 
DO 100 I=1.™ 
DH( I )=DERY( 1) *TT 
HCI) =HCI)+0H(T) 
100 CONTINUE 
T=T4+TT 
IF(teLTeII) GO TO 200 
WRITE( 60301) TeRO7) sH(7)2C(7) 
TICN)=T 
QS(N)=RC7) 
N=N+1 
IFC(II-GE.5) GO TO 188 
Il=I11+9 
GOor"Vo'*200 
188 II1=I11+10 
200 CONTINUE 
ND=K 
CALL DET3(TT2eQsUsND,TIER) 
CALL SE15(UsQesNDo IER) 
DO 700 1f=1.ND 
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IFCJIMeEQe2) GO TO 230 
UCT )=QC( I) 
GO TO 700 
230 VQ(I)=Q(1) 
700 CONTINUE 
KU=ND/5 
DO 124 TI=1+KU 
QU(1I)=U(1I*5) 
TUCI)=TACI*S5) 
WRITE( 62501) TUCI)»QU(I) 
124 CONTINUE 
301 FORMAT(C#®0'45X oF 5e153F 100-3) 
501 FORMATC#0*% 6 5XeF50els10XsFl00e3) 
CALL CGPL (TUse QUeTUsKU eo MK ol oleS3elsOeelOesl2ex0exles2es ALPHAs 6) 
350 CONTINUE 
MK=MK+41 
DO 300 LK=1 33 
READ(6235) CHsPPP 
35 FORMAT(F6e1sF12e9) 
MK=129 
DO 210 IV=1s2 
NN=CH+ND/ 10 
P=PPP 
DO 208 J=1.NN 
NI=10*J 
IN=1 
SUM=0 e 
IF(JeGTeND/10) NI=NOD 
IFC JeGTeCH) IN=10*J—-10*CH 
DO 201 JJ=IN.NI 
IFC IVeEQe2) GO TO 211 
SUM=SUM#U( JJ) 
GO TO 201 
211 SUM=SUM4UQ (JJ) 
201 CONTINUE 
Q( J) =0e1 *SUM*XP/SPRE 
208 CONTINUE 
TI(1)=le 
00 202 J=2.120 
TIC JV =TICJ—-1) 41. 
202 CONTINUE 
299 CONTINUE 
CALL CGPLCUT Ie Qe TI oNNoMK ol 010391 200210091 2020092059809 ALPHA,G) 
MK=MK+41 
210 CONTINUE 
CALL STATE 
CALL WAVE 
300 CONTINUE 
STOP 
END 
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SUBROUTINE STATE 


COMMON MKeSeALPHAsCHePPPelKs IVA 
REAL H(30) »DERY (30) »R( 30) 9 A(30) s RR( 30) ».BB( 30) C030) »MULTsDH(30) » 
1QS(€200) + TI( 200) sALPHA( 20) sN20FC 7) 0 QX( 200) sQY( 200) QZ( 200) 
M=7 
G=981le 
ALP=8.*G*S/V 
TT=e1l 
DO 60 I=1.M 
ACI )=€1-*91 «4/M 
60 CONTINUE 
K=1200 
FC=31.e 
DO 2 J=19M 
C(JJ=FC 
H( J)=06 
2 CONTINUE 
T=06 
IIl=1 
N=1 
DO 200 L=1>5K 
P=0. 
IF(TeLEeCH) P=PPP 
IF(PeEQe0-) GO TO 50 
DO SO I=1.5M 
R(T =ALP*61 eo * (HCI) *¥*36)/7C( 1) 
IFCR( I) -EQeO00.) GO TO SO 
QC=R(1)/761.6 
ABC=2 eS*ALOG( (100000. *P*V*H(I))7(01264*QC*QC) ) 
C(1)=24-+ABC 
IF(CC(CI)-GTeFC) CCII=FC 
SO CONTINUE 
DO 51 I=15M 
IF(PeEQeOe) C(II=24. 
R (1) =ALP*61-¢*( HCI) **32.)/C(T) 
51 CONTINUE 
DERY(1)=P—-R(1)7AC1) 
DO 70 I[=2sM 
DERY( 1) =P4( RC I-1-RCI)IZACT) 
70 CONTINUE 
DO 100 I=1.M: 
DH( I )=DERVCI)*TT 
H( 1) =HC I) +DHC I) 
100 CONTINUE 
199 T=T+TT 
IF(LeLTeII) GO TO 200 
WRITE(62301) TeR(7)5HO7),C(7) 
TICN)=T 
QS(N)I=R(7) 


ae ton 
i” 

“Ave ueeeanas ele, 

LOE VHOs TJUM + (OE) 30 (OE IGG + LOE IAA (OEIAg (OED Ae (OED YASG, € 
(00S 50+ (OOS) YH (OOS IXD o ET J5054 (OSAMA, (9083 E44 


| i + ; 
i 

li) 
ra 


en 


pan asa 

oe dT 0d 203d at 

| Me f=I_ 288 

CEIONG SERRE TIHD HS LOS han 

Liane 20 (20403. 399 be 

Whe PONCE! 

( (20820885 19 NCC TIHRVERN 0.0000 17920IR8BeS 


“Retr DD (24. 798199) 

re Hoe Ui . 

pare, esmtnna (20r0a a9 
(DN CAERE CL HDH e #9 ee oe: 


ed wl 12 
Bane a paNC rages a¥Aae Le 

ee Wen eet 
pre whi oF 


Al 


209 


D7 


N=N41 
IF(IIeGEeS) GO TO 188 
ITI=1L1i1¢+9 
GO TO 200 
188 II=11+10 
200 CONTINUE 
ND=N-1 
CALL CGPL(TIsQSoTI eoNDOoMKolsls4e1ls0eel0esl2es0022e59 Se sALPHAs6) 
MK=MK+1 
400 CONTINUE 
301 FORMAT(C® *& ,F200e1%93F20 eS5sF 3003) 
RETURN 
END 
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330 
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D8 
SUBROUTINE WAVE 


CCMMON MK oS sALPHAsCHsPPP.LKsIVA 
REAL NsGS(200)sTI(200) sHH( 200 )s ALPHA( 20) 
F (CH) =CALP*( H*¥*N) /VO)4tN*XAL PHT AX( HEX (N—-1 0) D—L 
L=91.4 

B=61l. 

G=981le 

NEW=e0093 

AA=.-00005 

N=3e 

TOT=120. 

TO=CH 

VO=PPP 

ALP=G*S/(3-*e0093) 

tt=1 

I=1 

T=e1l 

TS=(L/CALP*(VCX* (N—16)))) **C1 e/N) 
WRITE(6%100) TS 
FORMAT(*0%,10X9*TS=*,F5e2) 
IF(TOeLEeTS) GO TO 30 
IF(TeGTeTO) GO TO 77 
IF(TeLEeTS) Q=ALP*( VO XN) *(T &XN) *B 
IF(TeGTeTS) Q=L*VO*B 
IF((10*1I)-GTettL) GO TO 80 
QS(1)=Q 

TICI)J=T 

I=I+¢1 

T=Ttel 

LL=LL+1 

GO TO 76 

HO=(Q/( BX ALP) )**(1./N) 
WRITE(6:84) HO 
FORMAT( *0® o*FO=* Fl 0e7) 
HO=VOx#TS 

WRITE(6.84) HO 

GO TO 79 

IF(TeGTeTO) GO TO 330 

Q=BxALPx( (VOT) **N) 
IF((10*1I)-eGTe-LL) GO TO 31 
Qs(I)=Q 

TICI)=T 

I=1I¢1 

T=T+e1 

LL=LL+1 

GO TO 30 

HO=(Q/(B* ALP) )**(12/N) 

WRITE( 6484) HO 

HO=VO*«TO 

WRITE(€6e84) HO 
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101 
78 


81 


79 


99 


l2 


10 


13 


83 
82 


92 
90 
50 


TOS=L/ (ALP*( HOK*(N-10))) 
TP=(TOS-TO)/N 

WRITE(6e101) TP 

FORMAT °0%s10Xe* TP=* gF5e2) 
IF(TeGTe(TO+TP)) GO TG 79 
Q=Q 

IF((C10*I).GTetLtL) GO TO 81 
QS(1)=Q 

TICI)=T 

I=I+1 

T=T+el 

tL=LL+1 

GO TO 78 

TT=0e1 

T=T-.1 

IM=1 

CONTINUE 

TA=TT 

IF (TOcecLTeS02) TA=TPHTT 

J=0 

HO=HO 

H1=0.-01 

H2=( (HO*FCH1 ) )—CH1*FCHO)) ) (CF CHI )-FCHO) ) 
Z=ABS (CH2—-H1 )/ABSCH2 ) 
IF(ZeLEeAA) GC TO 10 
IF(JeEQe100) GO TO 11 
J=J+1 

HO=H1 

H1i=H2 

GO TO lt2 

WRITE(6 04) 

FORMAT(® *,232Xe*FAILED TO CONVERGE *) 
IFC IMeEQe2) GC TO 13 
IF(€C10*1I)-GTetLtL) GO TC 83 
Q=ALP*B*(H2**N) 

HH(I)=He 

QS(I)=Q 

TICI)=T+4TT 

T=I+1 

IM=2 

TT=TT+1. 

GO TQ 82 

TT=TT+el 

LL=LL+41 

IF( TT LEe(TOT-T)) GO TO 99 
NO=I 

DO 92 J=1sND 

WRITE (6090) TICJ) sQS(J) sHH( J) 
CONTINUE 

FORMATC® © 65X 0F10¢295Xe2F8 24) 
CONTINUE 
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IVA=IVA+1 


IFC IVAecEQe¢4) MK=—MK 
CALL CGPL(TIsQSeTIeNDeMK ol 91920120 eel0esl2ex0e92e528e0 sALPHAs6) 


RETURN 
END 
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MEASURED AND SIMULATED HYDROGRAPHS 
BY THE STATE VARIABLE MODEL 


RUNS NO. 1 TO 6, AND 10 TO 24 
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F-1 Confidence Intervals of Cumulative Hydrographs 


Cumulative outflows recorded during experiments performed on a 
laboratory catchment (single pulse rainfall inputs) were analysed 
statistically as described in Chapter VI. As a sample of the analyses, 
Table F-1 shows a computation of the 95% confidence intervals of the 
mean x (mean of 7 tests) of the accumulated outflow for the run number 
7. The confidence intervals were computed for Student's t = 2.447 and 


n = 7 according to the formula 


Kteskiv== (F-1) 
where s is the estimate of the standard deviation. The largest value 
of « was found to be + 5.49 cc which is less than the estimated maximum 
error of 6 cc given in Table 6-1 in Chapter VI. 

Table F-1 also contains values of accumulated outflow as predicted 
by the state variable model. It can be seen that the predicted values 
are within the confidence interval of the experimental cumulative hydro- 
graph. Because of the very good fit the cumulative hydrograph predicted 
by the state variable model was taken as a theoretical smooth curve 


representing the experimental data. 


F-2 Confidence Intervals of Instantaneous Outflow Hydrographs 


The instantaneous outflow rates were obtained from the measured 
cumulative outflow hydrographs by numerical differentiation. The 
computation of confidence intervals of these determined instantaneous 
outflow rates for the run number 7 is shown in Table F-2. 

Because numerical differentiation involves trunkation and round- 
off errors, it was necessary to establish the effect of differentiation 


on the confidence intervals given in Table F-2. 
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For this purpose, the theoretical curve of the cumulative outflow 
(predicted by the state variable model) was differentiated and the 
truncation error was determined (runoff error was neglected). The 
theoretical and experimental instantaneous outflow hydrographs were than 
compared to check the validity of the confidence intervals based on the 


differentiated experimental data. 


F-2-1 Determination of the Truncation Error 

Numerical differentiation was carried out on the IBM model 360-67 
computer using the subroutine DET3 described in the System/360 Scientific 
Subroutine Package (PL/I). The subroutine computes a vector Z = (2], 
Z25 see, Zn) Of derivative values, Gives SAVECtOrRYe=A(y1fay.2sennk,. yn) 
of function values whose components yj correspond to n equidistantly 
spaced argument values xj with xy - xy.] = h for i=l, 2, +. 25 0ns ihe 
resultant value zj involves a truncation error Ty where (Hildebrand, 1956) 


h2 
“a 96?) (G7) sce roleeex sd ie A= 1 


h2 
epece oy (3) Sei) eroixie rs aye; ) ie Se nel (F-2) 


I 
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” 
Boy) (En), Enelxy_p, X_] if 


Typical values of the theoretical curve of cumulative outflow, 
known at discreet time intervals with h = 0.1 sec, are shown in Table 
F-3. The third derivative, y (3), of the curve was determined from a 
relationship between differences and derivatives of a function (Conte, 
1965). In general 

AK f(x) = nk £€k) (¢) (F-3) 
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where AK £(x) is the kth difference of a function and £(k)(E) is the kth 


derivative of a function. The truncation error Ty in terms of differences 


is then given by 


Stab , €) elxy, x3] if i= 1 
hoes 
Ty =< - Ay tei) sea) sley-1, x771] 1f i= 2, ..., n-1 (F-4) 
3 
Aya) eh clx, 8 xe en 


3h 


Table F-3 also contains the computed first, second and third 
differences and the truncation errors given by Equation (F-4). The 
maximum error was found to be 0.0007 cc/sec, which is insignificant in 
comparison with the confidence limits given in Table F-2 and which 
ranged from 0.1 to 0.5 cc/sec. Since the theoretical curve of the 
instantaneous outflow lies within the confidence band of the experimental 


curve, the confidence band was accepted as valid. 
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SAMPLE CALCULATION OF THE TRUNCATION ERROR OF THE THEORETICAL 


Time 
(sec) 


15.0 


INSTANTANEOUS OUTFLOW HYDROGRAPH, RUN NO. 7 


Cumulative 
Discharge 


(ec) 


5.213 
5.354 
5.498 
5.644 
5.794 
5.946 


6.101. 


6.259 
6.421 
6.585 
6.753 
6.925 
TielOs 
7.280 
7.463 
7.650 
7.841 
8.036 
8.234 
8.438 


8.645 


Ay 


0.14082 
0.14365 
0.14651 
0.14941 
0.15234 
Q15532 
O51 5832 
0.16137 
0.16450 
0.16805 
0.17171 
0.17544 
Oe L7922 
0.18306 
0.18695 
0.19089 
0.19489 
0.19894 
0.20304 


0.20720 


TABLE F-3 


A2y 


0.00283 
0.00286 
0.00290 
0.00294 
0.00297 


0.00301 


0.00305 
000313 
0.00355 
0.00366 


0.00373 


0.00378 


0.00384 
0.00389 
0.00394 
0.00400 
0.00405 
0.00410 


0.00416 


A3y 


0.00004 
0.00004 
0.00004 
0.00004 
0.00004 
0.00004 
0.00009 
0.00042 
0.00011 
0.00006 
0.00006 
0.00006 
0.00005 
0.00005 
0.00005 
0.00005 
0.00005 


0.00005 


Truncation 


Error 


(cc/sec) 


A3y 
6h 


-0.00006 
-0.00006 
-0.00006 
-0.00006 
-0.00006 
-0.00006 
-0.00014 
-0.00070 
-0.00019 
-0.00010 
-0.00009 
-0.00009 
-0.00009 
-0.00009 
=0.00009 
-0.00009 
-0.00009 


-0.00009 
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APPENDIX G 
EXAMPLES OF THE PHASE PLANE REPRESENTATION 


OF SECOND ORDER SYSTEMS 
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G-l1 Determination of Phase Ipajectories by Analytical Solution for a 
Linear Conservative System. 


Consider a system described by the second-order differential 
equation 

24220 (6-1) 
Defining the state variables x, | 2, and x5 = z, Equation (G-1) may be 


written as 


X] = x9 os 
C= 
X2 = “Hy 
and 
an ks 
xy 7" x9 (G-3) 


In this case the variables in Equation (G-3) can be separated, thus 


XodXo = -K, dx] (G-4) 
and by integration the equation of trajectories is 
2 2 
Xo * x= C (G-5) 


Equation (G-5) defines phase trajectories which are concentric circles 
with radii r = c%, Some of the trajectories are shown in Figure G-l. 
The constant C is given by the initial conditions at time ty, x ,(to) and 
Ko (ta) - 

The system given by Equation (G-1) in the above example belongs 
to the class of second-order systems of the form zZ +.f£(z) = 0, which is 
known as that of conservative systems. It is so named because 2, which 
generally represents a measure of damping or energy dissipation in the 
system, is missing (Hsu, 1968). For this class of systems therefore, 
energy is conserved. This fact is illustrated by closed trajectories, 


such as the circles in Figure G-l, or can be seen by writing Equation 
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TRAJECTORIES 


Figure G-l. Phase trajectories for a system 


described by 2 +z = 0. 
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G 
(G-3) more generally as : 
dx2 _ -f(x1) (G-6) 
dx, x9 
Integrating Equation (G-6) yields 
4s. 
2 
*2 
sap Tate 8 £(x)dx = E (G-7) 
Xlo 


where E is the constant of integration. The first term of Equation 
(G-7) is seen to be a measure of the kinetic energy of the system, while 
the second term is a measure of the potential energy. xj, is the 
equilibrium state in x, with reference to which eneeseteneial energy is 
measured. The value of the constant E is thus a measure of the total 
energy; it is seen to be dependent on the initial condition. For each 


given value of E, Equation (G-7) gives the trajectory for the system. 


G-2 Determination of Phase Trajectories by the Method of Isoclines 


If the system differential equation cannot be solved easily then 
the trajectory must be found by numerical or graphical methods. One of 
the graphical methods is the method of isoclines. To describe the 


method in general terms, consider a system of the form 


xe = £7(x],*k2) 
(G-8) 


Xo = £5 (x), Xo) 


or 


Se ee a G-9 
x, 9X] PAG 4 ge od, oi 


(8-0) 
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There are three variables in Equation (G-9), X}> X9, and dx5/dx, = N(x), Xo). 
As defined, N(x], x») gives the slope of the local tangent of the trajectory 
passing through the point x), x9. If a specific numerical value is chosen 
for N, such that N = Ny, then 


EO(xyis x9) 


i. £1(x,, x) (G-10) 


Equation (G-10) is a function of x, and x9 only, thus defining a curve 

on the xj] - x9 plane, and this curve has the special property of being 
the locus of all points for which the slope of the phase trajectory is 
Nj}. The line is therefore called an isocline (Thaler, 1962). A suitable 
field of local tangents Ny, i=1, 2, ..., m, is then constructed 
assuming that the tangent is constant over a small region near the 

point. From this, the system trajectories can be constructed, if the 
initial conditions are given. The construction procedure is shown in 


Figure G-2. The following example serves to illustrate the method. 


G-2-2 Linear Non-conservative System 


Consider a system represented by the linear differential equation 
z+z2+22=0 (G-11) 
The state variables are again xj = z, and x9 = z. Equation (G-11) is 


then rewritten as 


x] = Xx 
: (G-12) 
Ray Ky 
or 
x =x] - x 
Lda int Conall (G-13) 
xX] pe 


from which 
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Figure G-2. 


INITIAL 
POINT A 


Construction of a trajectory (ABCD) using isoclines. 
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(G-14) 


Equation (G-14) defines the isoclines to be a family of straight lines 
passing through the origin of the phase plane. Figure G-3 shows isoclines 


and one particular trajectory for the system given by Equation (G-11). 


G-2-3 Nonlinear System. For nonlinear systems the isoclines are curves 
which usually require the computation of a number of points for each 


isocline. Consider the nonlinear system described by the equation 


2+ 22 "+ z= 0 | (G-15) 
from which 
Se OS (G-16) 
x] sae 
The isoclines are defined by 
x2 = ee (6-17) 


Family of isoclines given by Equation (G-17) and several phase trajectories 


of the system defined by Equation (G-15) are shown in Figure G-4. 
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Figure G-3. 


Phase plane diagram showing isoclines and a typical 
trajectory for a system described by xX + x +x 
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Figure G-4. Phase plane diagram for a system described 
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